笔记：Inductive Robust Principal Component Analysis

最新推荐文章于 2023-08-03 15:50:31 发布

XueShengke

最新推荐文章于 2023-08-03 15:50:31 发布

阅读量1.6k

点赞数

分类专栏：图像处理稀疏表示文章标签： inductive RPCA low-rank

图像处理同时被 2 个专栏收录

13 篇文章 8 订阅

订阅专栏

稀疏表示

12 篇文章 2 订阅

订阅专栏

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 \in R d, (1)

$\begin{equation} \tag{1} x = y + e, \quad x \in \mathbb{R}^d, \end{equation}$
其中

y $y$ 是低秩的子空间部分，

e $e$ 是误差项。给定任意的数据向量

x $x$ ，求解的目标是从原始数据中分离出向量

y $y$ 和误差

e $e$ 。

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

min U | | X - U U T X | | 2 F, s . t . U T U = I r, (2)

$\begin{equation} \tag{2} \min_U \ || X-UU^TX ||_F^2 \ , \ \ \mathrm{s.t.} \ U^T U = I_r \ , \end{equation}$
其中

Ir $I_r$ 是一个

r×r $r \times r$ 的单位矩阵，

||⋅||2F $||\cdot||_F^2$ 表示 Frobenious 范数，

r $r$ 是

U $U$ 的秩（

U $U$ 的列数）。这个优化问题可以有效地用奇异值分解（SVD）求出。假设

U∗ $U^*$ 是从训练数据

X $X$ 求出的解，有了一个数据

x $x$ ，它的主成分可以估计出

y = U * (U *) T x . (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 = [x_1, x_2, \cdots, x_n]$ ，RPCA 找到其主成分 $Y = [y_1, y_2, \cdots, y_n]$ ，通过如下的凸优化问题

min Y, 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∗)T $Y^* = U^* \Sigma^* (V^*)^T$ ；然后用

y=U∗(U∗)Tx $y = 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) 模型

min Z, E | | Z | | * + λ | | E | | ℓ, s . t . X = A Z + E, (5)

$\begin{equation} \tag{5} \min_{Z,E} \ ||Z||_* + \lambda ||E||_\ell, \ \ \mathrm{s.t.} \ X = AZ + E, \end{equation}$
其中

A $A$ 是数据空间的字典，

||⋅||ℓ $||\cdot||_\ell$ 根据具体情况选择合适的范数。如果

A=I $A = I$ ，那么 LRR 就变成了 RPCA 。子空间分割可以用 LRR 解决

min Z, E | | Z | * + λ | | E | | 2, 1, s . t . X = X Z + E, (6)

$\begin{equation} \tag{6} \min_{Z,E} \ ||Z|_* + \lambda ||E||_{2,1}, \ \ \mathrm{s.t.} \ X = XZ + E, \end{equation}$
其中区别是

A=X $A = 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 是可行的。首先，即使损坏是乱序的，也能存在一个线性投影 $P_0$ 把数据投影到子空间中，能正确地恢复出数据（即使不是完全精确的恢复）；
这里写图片描述
其次，两个高维的向量，通常是独立的，近似是相互正交 4，也就是说，损坏通常不在正确的子空间中。这样的情况下， $P_0$ 严格地从数据中去除损坏。只要有数据 $x$ ，其主成分就可以通过 $y = P_0 x$ 获得。

这里写图片描述

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

min P, E r a n k (P) + λ | | E | | 0, s . t . X = P X + 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$ 代替。重新构建的凸的优化问题

min P, E | | P | | * + λ | | E | | 1, s . t . X = P X + E . (8)

$\begin{equation} \tag{8} \min_{P,E} \ ||P||_* + \lambda ||E||_1, \ \ \mathrm{s.t.} \ X = PX + E \ . \end{equation}$
假设最优解

P∗ $P^*$ 已经得到。对于一个新的数据

x $x$ ，我们估计其主成分

y=P∗x $y = P^* x$ 和误差

e=x−P∗x $e = x - P^* x$ ，计算很快速、简便。

3.2 优化求解

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

min P, E | | P T | | * + λ | | E T | | 1, s . t . X T = P T X T + E T . (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^*$ 总是在

X $X$ 的各列分布的字空间中。

P∗ $P^*$ 可以被分解为

P∗=L∗(Q∗)T $P^* = L^* (Q^*)^T$ ，其中

Q∗ $Q^*$ 可以由正交化

X $X$ 的各列得到。于是，问题被转化为如下的形式

min L, E | | J | | * + λ | | E | | 1, s . t . X = L A + 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∗)TX $A = (Q^*)^T X$ ，

||⋅||∗ $||\cdot||_*$ 表示核范数，

||⋅||1 $||\cdot||_1$ 表示

ℓ1 $\ell_1$ 范数，

λ $\lambda$ 是权衡参数。

构建 Lagrangian 函数

L (L, E, J) = | | J | | * + λ | | E | | 1 + t r (Y T 1 (X - L A - E)) + t r (Y T 2 (L - J)) + μ 2 (| | X - L A - E | | 2 F + | | L - J | | 2 F), (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×n $Y_1,Y_2 \in \mathbb{R}^{m \times n}$ 分别是 Lagrange 乘子矩阵，

tr(⋅) $\mathrm{tr}(\cdot)$ 是迹函数，

μ $\mu$ 是惩罚项系数（非变量），

||⋅||F $||\cdot||_F$ 表示 Frobenius 范数。显然，优化目标函数可以给出

(L k + 1, E k + 1, J k + 1) Y 1 k + 1 Y 2 k + 1 μ = arg min L, E, J L (L, E, J, Y 1 k, Y 2 k), = Y 1 k + μ (X - L A - E), = Y 2 k + μ (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），

L k + 1 E k + 1 J k + 1 = arg min L L (L, E k, J k, Y 1 k, Y 2 k), = arg min E L (L k + 1, E, J k, Y 1 k, Y 2 k), = arg min J L (L k + 1, E k + 1, J, Y 1 k, Y 2 k) . (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) = s i g n (x) \cdot 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，

L k + 1 = ((X - E k) A T + J k + Y 1 k A T - Y 2 k μ) (I + A A T) - 1; E k + 1 = S λ μ (X - L k + 1 A + Y 1 k μ); (U, Σ, V) = s v d (L k + 1 + Y 2 k μ), J k + 1 = U S 1 μ (Σ) V T; (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+1 $L_{k+1},E_{k+1},J_{k+1}$ 的过程中，并没有考虑这两项，而是用减去

Ek $E_k$ 的较小的元素值

Sλ/μ(Ek) $S_{\lambda/\mu}(E_k)$ ，和减去

Jk $J_k$ 的部分（较小的）奇异值

S1/μ(Σ) $S_{{1}/{\mu}}(Σ)$ 来代替，直接求闭式解表达式。
优化步骤的迭代停止条件：

| | X - L k + 1 A - E k + 1 | | \infty < ε a n d | | Z - J k + 1 | | \infty < ε . (20)

$\begin{equation} \tag{20} ||X-L_{k+1}A-E_{k+1}||_\infty < \varepsilon \ \mathrm{and} \ ||Z-J_{k+1}||_\infty < \varepsilon. \end{equation}$
其中

||⋅||∞ $||\cdot||_\infty$ 表示最大范数，定义为矩阵中最大的元素绝对值

||A||∞=maxi,j{|aij|} $||A||_\infty = \max_{i,j} \left\{ |a_{ij}| \right\}$ 。

4 实验

略

Appendix

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

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

$\partial L \partial L \partial L \partial L L k + 1 = \partial \partial L {t r (Y T 1 k (X - L A - E k)) + t r (Y T 2 k (L - J k)) + μ 2 (| | X - L A - E k | | 2 F + | | L - J k | | 2 F)} = - Y 1 k A T + Y 2 k - μ (X - L A - E k) A T + μ (L - J k) = μ L (I + A A T) + (Y 2 k - Y 1 k A T) + μ ((E k - X) A T - J k), = 0, = ((X - E k) A T + J k + Y 1 k A T - Y 2 k μ) (I + A A T) - 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}$
求解 $E_{k+1}$ , 当 $\frac{\partial \mathcal{L}}{\partial E} = 0$ 时，（ $||J||_*$ 与 $E$ 无关， $\lambda ||E||_1$ 用 $S_{\lambda / \mu} (\cdot)$ 近似）

$\partial L \partial E \partial L \partial E E k + 1 = \partial L \partial E {t r (Y T 1 k (X - L k + 1 A - E)) + μ 2 | | X - L k + 1 A - E | | 2 F} = - Y 1 k + μ (E + L k + 1 A - X), = 0, = S λ μ (X - L k + 1 A + Y 1 k μ) .$ $\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}$
求解 $J_{k+1}$ , 当 $\frac{\partial \mathcal{L}}{\partial J} = 0$ 时，（ $\lambda ||E||_1$ 与 $J$ 无关， $||J||_*$ 用 $S_{1/ \mu} (\Sigma),\ \Sigma = \mathrm{svd}(\cdot)$ 近似）

$\partial L \partial J \partial L \partial J (U, J k + 1 = \partial L \partial J {t r (Y T 2 k (L k + 1 - J)) + μ 2 | | L k + 1 - J | | 2 F} = - Y 2 k + μ (J - L k + 1), = 0, Σ, V) = s v d (L k + 1 + Y 2 k μ), = U S 1 μ (Σ) V T .$ $\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 ↩