RPCA模型及GoDec算法

一、基础知识

1. 在线性代数中，秩的物理意义是矩阵中最大的不相关的向量个数，因此秩可以用来度量相关性；而在图像处理中，秩就可以理解为图像中所包含的信息的丰富程度。图片中比较突兀的成分，会增加图像矩阵的秩，一幅清晰的自然图像其数据矩阵往往是低秩或者近似低秩的，这是因为其中的图像信息有很大的相关性，但如果图像中引入噪声，那么存在随机幅值任意大但是分布稀疏的误差破坏了原有数据的低秩性。
2. 低秩矩阵：如果$A$是一个$m\times n$的矩阵，$rank(A)$是矩阵$A$的秩，如果$rank(A)$远小于$m$和$n$，则我们称$A$是低秩矩阵。低秩矩阵说明这个矩阵包含了大量的冗余信息，秩越低表示数据冗余性越大。
3. 低秩是指矩阵的秩较小
4. 稀疏是指矩阵中非零元素的个数少，即大部分元素都为0，只有小部分元素不为0，用$l_0$范数度量
5. 矩阵$L$的0范数：矩阵$L$中非零元素的个数；非零元素个数越少，$L$越稀疏
6. 矩阵$L$的核范数(nuclear norm)：矩阵$L$的奇异值之和；核范数越小，可近似认为矩阵$L$的秩越低
7. 矩阵$L$的1范数：矩阵$L$中各列元素绝对值的和的最大值；这个和很小，趋近于0的时候，可近似认为非零元素个数很少。
   $$||L|{|}_1=\underset{j}{\max }\sum _{i=1}^m|a_{ij}|$$
8. P问题(Polynomial Solvable)：可以在多项式时间内解决的问题；时间复杂度如$（n^2,n^4,n(log(n))）$都是多项式时间
9. NP问题：多项式复杂程度的非确定性问题
   NP难问题的进一步理解

---

二、鲁棒主成分分析（RPCA）

概述

低秩矩阵恢复模型（LRMR）主要有鲁棒主成分分析（RPCA）、矩阵补全（matrix completion，MC）和低秩表示（low-rank representation，LRR）等三类模式。

低秩矩阵恢复问题，也被称为低秩稀疏分解问题，又被称为鲁棒主成分分析问题，可描述为从一个已知矩阵 $\mathbf{D}\in {\mathbf{R}}^{\mathbf{m}\times \mathbf{n}}$恢复出一个低秩矩阵$\mathbf{L}\in {\mathbf{R}}^{\mathbf{m}\times \mathbf{n}}$和一个稀疏矩阵$\mathbf{S}\in {\mathbf{R}}^{\mathbf{m}\times \mathbf{n}}$。¹
$$D=L+S$$

数学模型

低秩矩阵恢复问题的数学模型可表示为
$$\underset{L,S}{\min }rank(L)+\lambda ||S|{|}_{0}s.t.D=L+S\text{\hspace{1em}(1)}$$
其中$rank(L)$表示矩阵的秩，$||S|{|}_0$表示矩阵$S$非零元的个数。
 
表示的含义是当低秩矩阵$L$的秩和稀疏矩阵$S$的0范数均取得最小值的时候，所得到的对应的$L$和$S$.

上式中秩函数、零范数都是非凸非线性的组合优化函数，模型1是非凸NP-难问题，直接求解很难。需要在一定条件下对其松弛，方可进行优化。

用矩阵的核函数来近似矩阵的秩，用矩阵的1范数来近似0范数³。由范数知识可知，核函数是秩函数的凸包，1范数是0范数的凸包，所以上述NP-难问题松弛后可转化为凸优化问题：
$$\underset{L,S}{\min }||L|{|}_{\ast }+\lambda ||S|{|}_{1}s.t.D=L+S\text{\hspace{1em}(2)}$$
式（2）中的核范数是指矩阵$D$的奇异值的和，主要用来约束$D$的低秩属性，$L_1$范数用来约束误差图像矩阵$S$的稀疏性。²

基于凸优化松弛模型 2，已经有很多算法被提出，其中包括増广拉格朗日乘子法（augmented Lagrange multiplier，ALM）、加速的近似梯度法（accelerated proximal gradient，APG）和交替方向法（alternating direction method of mutipliers，ADMM）。¹

矩阵低秩分解理论
矩阵分解在推荐系统中的应用

---

三、GoDec算法⁴

GoDec原理

GoDec算法是低秩稀疏领域中一种新型的快速近似分解方法。基本思想是将一个原始矩阵快速分解成低秩矩阵、稀疏矩阵和噪声矩阵三部分。且可以指定$X-S$到$L$的低阶近似和$X-L$到$S$的稀疏近似。利用预先设定的低秩参数$r$和稀疏参数$k$，可以将不具备严格低秩稀疏结构的矩阵进行低秩稀疏近似分解。

其数学模型如下：
$$D=L+S+N，rank(L)\le r，card(S)\le k\text{\hspace{1em}(3)}$$
式中，$rank(L)$表示对矩阵$L$求秩，$card(S)$表示矩阵$S$的稀疏度，$k$值越小，矩阵$S$越稀疏。

求解过程

对于上式中的分解问题，可以通过最小化分解误差来解决，即
$$\{\begin{array}{l}\underset{L,S}{\min }||D-L-S|{|}_F^2\\ s.t.rank(L)\le r，card(S)\le k\end{array}\text{\hspace{1em}(4)}$$
式中，$||\cdot |{|}_F$表示矩阵的$F$范数，即使低秩稀疏分解后的重构误差最小。

式(4)优化问题可以分解成下面两个子问题进行交替求解，直至收敛。
$$\{\begin{array}{l}{L}_t=arg\underset{rank(L)\le r}{\min }||D-L-S_{t-1}|{|}_F^2\\ \\ S_t=arg\underset{card(S)\le k}{\min }||D-L_t-S|{|}_F^2\end{array}\text{\hspace{1em}(5)}$$
虽然这两个子问题都是非凸约束问题，但它们存在全局解$L_t$和$S_t$。由上式可以看出，$L_t$和$S_t$可以通过互相迭代求解，类似于围棋中的双方博弈，这也是算法名称中“Go”的由来。
 
式（5）中子问题中的$L_t$和$S_t$都可以通过设置硬阈值的方式求得，具体如下：
$$\{\begin{array}{l}{L}_t={\sum }_{i=1}^r{\lambda }_i{U}_i{V}_i^T，svd(X-S_{t-1})=U\Lambda V^T；\\ \\ S_t=P_{\Omega }(X-L_t),\Omega :|(X-L_t{)}_{i,j\in \Omega }|\ne 0\\ and\ge |(X-L_t{)}_{i,j\in \bar{\Omega }}|,|\Omega |\le k.\end{array}\text{\hspace{1em}(6)}$$
$L_t$通过$X-S_{t-1}$的奇异值硬阈值（前r个奇异向量）更新；$S_t$通过$X-L_t$的entry-wise 硬阈值（数值从大到小排列前k个元素组成的非零子集）更新.
 
由于SVD通常很耗时，需要$min(m{n}^2,m^{2}n)$次计算，所以考虑用双边随机投影（Bilateral Random Projection，BRP）来计算式(5)，相比于传统计算方法SVD，运算速度明显提高，且很大程度地将计算复杂度降低。

为了快速分解出低秩矩阵和稀疏矩阵，根据给定的秩$r$，利用随机矩阵$A_1\in R^{n\times r}$和$A_2\in R^{m\times r}$构造密矩阵$D$的双边随机投影（bilateral random projection，BRP）$Y_1$和$Y_2$，即
$$Y_1=D{A}_1，Y_2=D^T{A}_2\text{\hspace{1em}(6)}$$
式中，$D$为$m\times n$(m≥n)的矩阵。
 
可以快速获得秩为$r$的低秩矩阵$L$
$$L=Y_1(A_2^T{Y}_1{)}^{-1}{Y}_2^T\text{\hspace{1em}(7)}$$
计算$L$包括一个$r\times r$矩阵的求逆和三个矩阵的相乘。对于密矩阵$D$，可以得到BRP的时间复杂度为2mnr，则得到矩阵$L$的必要时间复杂度为$r^2(2n+r)+mnr$，计算量远低于使用SVD分解，效率更高。

当矩阵$D$的奇异值衰减缓慢时，利用模型优化（Power Scheme）对矩阵进行优化。使用$\bar{D}=(D{D}^T{)}^{q}D$来代替$D$，$\bar{D}$的奇异值衰减速度要比$D$快，因此计算$\bar{D}$的BRP。特别地，${\lambda }_i(\bar{D})={\lambda }_i(\bar{D}{)}^{2q+1}$（${\lambda }_i$为矩阵的第$i$个最大奇异值），$D$和$\bar{D}$有相同的奇异向量。$\bar{D}$的BRP为
$$Y_1=\bar{D}{A}_1，Y_2={\bar{D}}^T{A}_2\text{\hspace{1em}(8)}$$
根据式(7)可以得到优化后的低秩矩阵$\bar{L}$
$$\bar{L}=Y_1(A_2^T{Y}_1{)}^{-1}{Y}_2^T\text{\hspace{1em}(9)}$$
为了得到矩阵$D$秩为$r$的低秩近似，对$Y_1$和$Y_2$进行$QR$分解
$$Y_1=Q_1{R}_1，Y_2=Q_2{R}_2\text{\hspace{1em}(10)}$$
最终得到矩阵$D$的低秩近似$L$为
$$L=(\bar{L}{)}^{\frac{1}{2q+1}}=Q_1[R_1(A_2^T{Y}_1{)}^{-1}{R}_2^T{]}^{\frac{1}{2q+1}}{Q}_2^T\text{\hspace{1em}(11)}$$
式中，$q$代表幂指数，可以通过增加$q$值来减少式(7)的误差。
 
为了精确恢复低秩矩阵$L$，使用右随机投影$Y_1$构建左投影矩阵$A_2$，并计算左随机投影$Y_2$。当通过$Y_1=D{A}_1$得到$Y_1$后，更新$A_2=Y_1$，计算左随机投影$Y_2=D^T{Y}_1$，然后更新$A_1=Y_2$，再计算右随机投影$Y_1=D{Y}_2$。通过不断更新式(7)中的$Y_1$和$Y_2$，得到更精确的低秩矩阵$L$。

然后可以得到稀疏矩阵$S$为
$$S=P_{\Omega }(D-L)\text{\hspace{1em}(12)}$$
式中，$P_{\Omega }(\cdot )$表示元素集$\Omega$对矩阵的采样投影，$\Omega$为$|D-L|$的前$k$个最大元素非零子集。

改进初始化矩阵

原始GoDec算法在进行分解时，需要初始化低秩矩阵$L_0$为原始矩阵$D$，稀疏矩阵$S_0$为零矩阵，表示成
$$\{\begin{array}{l}{L}_0=D\\ S_0=0\end{array}\text{\hspace{1em}(13)}$$
为了进一步提高计算速度，减少迭代次数，这里使用原矩阵每一行的中值替换该行所有值来逼近低秩矩阵$L_0$，用原矩阵减去低秩矩阵后剩余的矩阵作为稀疏矩阵$S_0$的初始值，则改进后的初始化矩阵为
$$\{\begin{array}{ll}{L}_0=med(D_i){1\text{1}1}_{1\times n}，& i=1,2,\cdots ,m\\ S_0=D-L_0& \end{array}\text{\hspace{1em}(14)}$$
式中，$med(D_i)$表示矩阵$D$第$i$行的中值，${1\text{1}1}_{1\times n}$表示$1$行$n$列的单位向量。

与RPCA的区别

GoDec（Go Decomposition）是Tao等人在RPCA的基础上提出的一种矩阵分解算法。与RPCA不同，它不再要求原始数据矩阵$X$中的噪声矩阵满足稀疏条件，而是对X进行了更精细的划分，从而实现矩阵分解。从数学表达式上讲，它将原始数据$X$分成了低秩矩阵$L$、稀疏矩阵$S$和噪声矩阵$G$，即$X=L+S+G$。此外，GoDec与RPCA还有以下三点不同：

1. 由式（2）可以看出，RPCA不需要提前设定低秩矩阵的秩，它通过最小化低秩矩阵的核范数与稀疏矩阵的$l_1$范数的和求得最优解，因此最终矩阵$L$的秩和矩阵$S$的稀疏度由凸优化算法自动求得。而DoDec可以在求解之前提前设定参数$rank(L)$和$card(S)$来约束矩阵$L$的秩以及$S$的稀疏程度，从而控制模型复杂度。
2. 由于GoDec可以控制$S$的支撑集，因此它可以用来解决矩阵补全（Matrix Completion）问题。
3. 与RPCA相比，GoDec通过对低秩成分的随机近似计算显著提高了运算效率，大大减少了运行时间。

---

四、参考文献

[1]. 求解低秩矩阵恢复问题的一种非单调交替梯度方向法_闫喜红

[2]. 基于判别低秩矩阵恢复和协同表示的遮挡人脸识别_孙雨浩

[3]. 非凸的鲁棒主成分分析模型及其应用_王宇

[4]. 基于GoDec的探地雷达信号处理方法研究_郑文俭