引言

在过去的十多年间，数据分析行业发生了革命性的变化。目前，用户需要处理的数据在尺寸、大小和复杂性上都呈爆炸式增长趋势^[1]，所以如何通过识别低维结构从海量数据中提取有用的信息，已成为一个迫切需要解决的问题。

主成分分析(PCA)是一种用于发现数据中低维结构的通用技术，它假设数据是从一个单一的低维数据中提取的高维空间的仿射子空间。PCA作为最简单、最流行的降维工具，在计算机视觉与图像处理等许多领域得到广泛的应用^[2-4]。然而，PCA对于含野点和稀疏大噪声的数据非常敏感，为解决此问题，研究者们提出了鲁棒主成分分析(RPCA)方法^[5]。RPCA方法通过最小化核范数和1范数的组合问题，将观测矩阵分离为低秩部分和稠密小噪声部分。RPCA方法已广泛应用于图像去噪、视频处理、网页搜索和生物信息等领域^[6-8]，但当数据被两种以上混合噪声污染时，RPCA的表现不太理想，故本文提出一种可处理含有混合噪声的模型算法即广义鲁棒主成分分析(GRPCA)，通过最小化核范数、1范数和2, 1范数的组合问题，从观测矩阵中分离出低秩部分和混合噪声部分，采用随机排序的交替方向乘子法(RPADMM)对GRPCA求解，并选取机器学习常用网站UCI中的一个垃圾邮件数据集(spambase)^[9-10]来验证GRPCA算法的性能。

1 鲁棒主成分分析

广义鲁棒主成分分析是由RPCA衍生出的一种新型降维算法，RPCA由PCA发展而来，可以用于低秩矩阵恢复，即从一个被稀疏大噪声污染的低秩矩阵中恢复出低秩部分，并分离出噪声。

PCA可以通过约束优化模型(1)求出低秩矩阵A的最优估计

$ \left\{ \begin{array}{l} \mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{A}}, \mathit{\boldsymbol{E}}} {\left\| \mathit{\boldsymbol{E}} \right\|_F}\\ {\rm{s}}{\rm{.t}}{\rm{.rank}}\left( \mathit{\boldsymbol{A}} \right) \le r, \mathit{\boldsymbol{D}} = \mathit{\boldsymbol{A}} + \mathit{\boldsymbol{E}} \end{array} \right. $

(1)

式中，‖·‖_F为Frobenius范数。设X是一个n维方阵，则$ {\left\| \mathit{\boldsymbol{X}} \right\|_F} = {\left( {\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{{({a_{i, j}})}^2}} } } \right)^{\frac{1}{2}}}$为矩阵元素的平方和的算术平方根。对观测矩阵D进行奇异值分解(SVD)即可求出最优解。

当D受到稀疏大噪声的干扰时，PCA的假设条件不满足，无法正常工作。此时，恢复出低秩矩阵A就变成如式(2)的非凸优化问题

$ \left\{ \begin{array}{l} \mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{A}}, \mathit{\boldsymbol{E}}} {\rm{rank}}\left( \mathit{\boldsymbol{A}} \right) + \lambda {\left\| \mathit{\boldsymbol{E}} \right\|_0}\\ {\rm{s}}{\rm{.t}}.\mathit{\boldsymbol{D}} = \mathit{\boldsymbol{A}} + \mathit{\boldsymbol{E}} \end{array} \right. $

(2)

式中，‖·‖₀为0范数，表示矩阵中非0元素的个数。式(2)的优化问题经松弛变为一个新的凸优化问题^[5]，即RPCA

$ \left\{ \begin{array}{l} \mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{A}}, \mathit{\boldsymbol{E}}} {\left\| \mathit{\boldsymbol{A}} \right\|_*} + \lambda {\left\| \mathit{\boldsymbol{E}} \right\|_1}\\ {\rm{s}}{\rm{.t}}.\mathit{\boldsymbol{D}} = \mathit{\boldsymbol{A}} + \mathit{\boldsymbol{E}} \end{array} \right. $

(3)

式中‖·‖₁为1范数，则$ {\left\| \mathit{\boldsymbol{X}} \right\|_1} = \left( {\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {|{x_{i, j}}|} } } \right) $表示矩阵各元素绝对值之和。式(3)的RPCA模型可用迭代阈值算法(IT)、交替方向乘子法(ADMM)^[6]或近端梯度算法(APG)^[7]等算法求解。

由于污染矩阵的噪声可能不唯一，RPCA算法无法处理被两种不同噪声污染的情况，所以本文在RPCA算法基础上加入新的噪声矩阵H，并采用RPADMM对GRPCA求解，该方法的收敛性已被证明^[11]。

2 广义鲁棒主成分分析(GRPCA)及其算法 2.1 模型及求解

当D同时被稀疏大噪声(如椒盐噪声)和稠密小噪声(如高斯噪声)干扰时，RPCA模型无法识别稠密小噪声, 因此本文在模型(3)中引入了第三项l_{2, 1}范数来表示稠密的小噪声。l_{2, 1}范数定义如下。

设X是一个n维方阵，且$ {\left\| \mathit{\boldsymbol{X}} \right\|_{2, 1}} = \sum\limits_{i = 1}^n {} \sqrt {\sum\limits_{j = 1}^n {x_{_{i, j}}^{^2}} } $, 此时RPCA模型可改写为

$ \left\{ \begin{array}{l} \mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{A}}, \mathit{\boldsymbol{E}}} {\left\| \mathit{\boldsymbol{A}} \right\|_*} + \lambda {\left\| \mathit{\boldsymbol{E}} \right\|_1} + \gamma {\left\| \mathit{\boldsymbol{H}} \right\|_{2, 1}}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\mathit{\boldsymbol{D}} = \mathit{\boldsymbol{A}} + \mathit{\boldsymbol{E}} + \mathit{\boldsymbol{H}} \end{array} \right. $

(4)

模型(4)的增广拉格朗日函数为

$ \left\{ \begin{array}{l} L\left( {\mathit{\boldsymbol{A}}, \mathit{\boldsymbol{E}}, \mathit{\boldsymbol{H}}, \mathit{\boldsymbol{Y}}, \mu } \right) = {\left\| \mathit{\boldsymbol{A}} \right\|_*} + \lambda {\left\| \mathit{\boldsymbol{E}} \right\|_1} + \\ \;\;\;\;\;\;\;\;\;\gamma {\left\| \mathit{\boldsymbol{H}} \right\|_{2, 1}} + < \mathit{\boldsymbol{Y}}, \mathit{\boldsymbol{D}} - \mathit{\boldsymbol{A}} - \mathit{\boldsymbol{E}} - \mathit{\boldsymbol{H}} > + \\ \;\;\;\;\;\;\;\;\;\mu \left\| {\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{A}} - \mathit{\boldsymbol{E}} - \mathit{\boldsymbol{H}}} \right\|_{_F}^{^2} \end{array} \right. $

(5)

由于式(5)含有3个变量(A，E，H)，ADMM求解不再收敛，故本文采用RPADMM对其求解。RPADMM与ADMM的区别在于，每一步迭代的顺序是不确定的，取决于迭代前产生的一组随机数即(1，2，3)的随机排序。随机数为(2，3，1)表示先求第2个变量，再求第3个变量，最后求第1个变量，即先固定A、H同时更新E，再固定A、E同时更新H，最后固定E、H同时更新A。

GRPCA求解迭代格式如下。

固定E、H，更新A时，有

$ {\mathit{\boldsymbol{A}}_{k + 1}} = {\mathit{\boldsymbol{D}}_{{\mu ^{ - 1}}}}(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{E}} - \mathit{\boldsymbol{H}} + {\mu ^{ - 1}}\mathit{\boldsymbol{Y}}) $

固定A、H，更新E时，有

$ \begin{array}{l} {\mathit{\boldsymbol{E}}_{k + 1}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{E}} \lambda {\left\| \mathit{\boldsymbol{E}} \right\|_1} + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{A}} - \mathit{\boldsymbol{E}} - \mathit{\boldsymbol{H}} + {\mu ^{ - 1}}\mathit{\boldsymbol{Y}}} \right\|_{_F}^{^2} = \\ {S_{\frac{\lambda }{\mu }}}(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{A}} - \mathit{\boldsymbol{H}} + {\mu ^{ - 1}}\mathit{\boldsymbol{Y}}) \end{array} $

固定A、E，更新H时，有

$ \begin{array}{l} {\mathit{\boldsymbol{H}}_{k + 1}} = \arg \mathop {\min }\limits_F \gamma {\left\| \mathit{\boldsymbol{H}} \right\|_{2, 1}} + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{A}} - \mathit{\boldsymbol{E}} - \mathit{\boldsymbol{H}} + {\mu ^{ - 1}}\mathit{\boldsymbol{Y}}} \right\|_{_F}^{^2} = \\ {\left( {\gamma \mathit{\boldsymbol{G}} + \mu \mathit{\boldsymbol{I}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{Y}} + \mu \left( {\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{A}} - \mathit{\boldsymbol{E}}} \right)} \right) \end{array} $

式中$ \mathit{\boldsymbol{G}} = {\rm{diag}}\left( {\frac{1}{{{{\left\| {{\mathit{\boldsymbol{H}}^i}} \right\|}_2}}}} \right) $是由H决定的对角矩阵。

GRPCA算法步骤如下。

1) 初始化参数Y₀=0，μ₀>0，ρ>0，γ=0.5，k=0，A₀=0，E₀=0，H₀=0

2) if 收敛条件满足，输出结果

3) else

    for k=0到最大迭代次数

      生成(1, 2, 3)的一个随机排序σ

4)    遍历σ

5)    if σ(i)=1

        A_k+1=D_μ^-1(D-E-H+μ^-1Y)

6)    else if σ(i)=2

        $ {\mathit{\boldsymbol{E}}_{k + 1}} = {S_{\frac{\lambda }{\mu }}}(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{A}} - \mathit{\boldsymbol{H}} + {\mu ^{ - 1}}\mathit{\boldsymbol{Y}})$

7)     else σ(i)=3

         H_k+1= (γG+μI)^-1(Y+μ(D-A-E))

8)       end if

9)    Y_k+1=Y_k+μ_k(D－A_k－E_k－H_k)

       μ_k+1=ρμ_k, k=k+1

10)  end for

    end if

11) 输出结果(A^*, E^*, H^*)

2.2 算法复杂度分析

在求解模型(4)的GRPCA算法中，运算量主要集中在奇异值分解和矩阵求逆中。假设输入矩阵的维度为R^m×n，每一次循环中，步骤5)和步骤6)奇异值分解的运算量为O((m+n)³)，步骤7)矩阵求逆的运算量为O((m+n)³)。此外GRPCA算法循环内的多次乘法加法和截断阈值的运算量为O(5n(m+n)+14n)。因此，GRPCA算法的总运算量为O((m+n)³+I(5n(m+n)+14n))，其中I为循环的次数。

3 实验验证及结果分析

所有实验运行环境均为Matlab R2016b版，机器配置为ThinkCentre M8400t-N000，内存4 GB，频率3.4 GHz，操作系统为Windows 7。

3.1 数据来源

本文数据集来自UCI中的Spambase Data Set。Spambase是一个4 601×58的矩阵，包含4 601个邮件，每一行代表一个邮件。其中前57列为邮件的属性，第58列为邮件的标签(1表示垃圾邮件0表示非垃圾邮件)。本文采用分类精确度(accuracy)来度量实验的有效性。

3.2 实验设置

实验分为两大部分。

(1) 将本文提出的GRPCA与kNN算法结合构成垃圾邮件分类的一个新型算法：先用GRPCA对邮件属性进行去冗余预处理，再用kNN分类。并与经典的降维算法PCA、RPCA和kNN结合所形成的垃圾邮件分类算法比较。

(2) 将GRPCA与SVM算法结合构成垃圾邮件分类的一个新型算法:先用GRPCA算法对邮件属性进行去冗余预处理，再用SVM算法分类。并与经典的降维算法PCA、RPCA和SVM结合所形成的垃圾邮件分类算法比较。实验流程如图 1。

3.3 结果分析

3.2节中两部分实验所得结果分别如表 1、图 2所示。从表 1可以看出，当固定kNN算法的k值时，经过GRPCA去冗余后再分类的精确度均达到了85.5%以上，比直接使用kNN的精确度要高出2.91%~4.79%；相对经典降维算法PCA和RPCA预处理后的结果的精确度分别高出2.91%~4.79%和1.22%~1.85%，并且3种方法所用时间几乎无差别。从图 2可以看出，在SVM算法中，经过GRPCA去冗余处理后的精确度高达89.52%，比直接使用SVM的精确度要高出2.76%，比PCA和RPCA预处理后的精确度分别高出2.76%和3.89%。由于3种方法所需时间均较长(大于6 h)，故不作时间对比。

此外，还可以看出，在PCA去冗余处理前后精确度几乎没有发生变化，这说明属性矩阵受到稀疏大噪声的干扰，使PCA无法正确降维，从而证明本文在邮件过滤时对属性矩阵进行去冗余处理是必要的。

4 结束语

本文提出了一种新的广义鲁棒主成分分析(GRPCA)算法，通过最小化核范数、1范数和2, 1范数的组合问题，可以分离出被混合噪声污染的低秩矩阵。选用垃圾邮件来检验算法，结果表明本文算法在提高垃圾邮件的精确度上是有效的。