稀疏低秩矩阵分解是近年来机器学习和数据分析等领域的热点问题，是高维数据处理的基础和核心问题之一。在鲁棒主成分分析^[1]、矩阵恢复^[2]和图像处理^[3]等领域均具有重要的应用价值。

假设D∈R^m×n是待分解矩阵，稀疏低秩矩阵分解的数学模型可描述为：

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{L,S \in {R^{m \times n}}} {\rm{rank}}\left( L \right) + \lambda {{\left\| S \right\|}_0}}\\ {{\rm{s}}.\;{\rm{t}}.\;\;\;\mathit{L} + S = D。} \end{array} $

(1)

其中: L为D的低秩成分；rank(L)表示L的秩；S为D的稀疏成分；用L₀模度量S的稀疏度，正则化参数λ>0用于平衡L的低秩程度和S的稀疏程度。

由于秩函数和L₀模的非凸不连续性，稀疏低秩矩阵分解的原始模型(1)的求解是NP-难问题。受压缩感知研究领域的启发，Candes等^{[1, 2, 4]}将模型(1)松弛到如下的凸模型，并且给出了两模型相应的等价性条件：

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{L,S \in {R^{m \times n}}} {{\left\| L \right\|}_ * } + \lambda {{\left\| S \right\|}_1}}\\ {{\rm{s}}.\;{\rm{t}}.\;\;\;\mathit{L} + S = D。} \end{array} $

(2)

其中：${{\left\| L \right\|}_{*}}$表示L的核范数，即所有奇异值之和；${{\left\| S \right\|}_{1}}$为众所周知的L₁模。近几年发展的非光滑凸优化算法如ADMM^[5-6]等均可以高效地求解替代模型(2)。

尽管模型(2)在稀疏低秩矩阵分解中已经取得了成功的应用，但其实际效果仍有巨大的提高空间。首先，模型(2)与模型(1)的等价性条件在实际问题中往往并不满足，因此模型(2)只是模型(1)的凸近似，二者之间的固有模型误差难以修复；其次，压缩感知^[7]、图像处理^[8]和矩阵恢复^[9]等领域的研究成果均已表明，采用非凸函数惩罚稀疏性优于目前广泛使用的L₁模。文献^[10]将一类非凸稀疏罚函数应用于压缩感知问题，发现其可进一步诱导向量的稀疏性，并且易于设计高效的数值算法。

本文将文献[10]中所采用的非凸稀疏罚函数推广到矩阵情形，分别用其惩罚低秩矩阵和稀疏矩阵的稀疏度，提出了稀疏低秩矩阵分解模型(1)的新的非凸近似模型。考虑到目标函数与约束条件关于L和S的可分性，本文基于ADMM算法框架提出了新模型的数值求解算法，并且给出了其迭代子问题最优解的具体形式，保证了所设计算法的高精度与低时间复杂度特性。数值实验结果表明，在稀疏度较小的情况下，本文的非凸模型相较于广泛使用的L₁模型(2)具有一定的优势。

1 非凸近似模型 1.1 阈值算子与非凸罚函数

对于矩阵函数f:R^m×n→R和参数β>0，其迫近算子prox_βf:R^m×n→R^m×n定义为：

$ {\rm{pro}}{{\rm{x}}_{\beta f}}\left( X \right) = \mathop {\arg \min }\limits_{Y \in {R^{m \times n}}} \left\{ {\frac{1}{2}\left\| {Y - X} \right\|_F^2 + \beta f\left( Y \right)} \right\}。$

(3)

稀疏低秩矩阵分解的凸近似模型(2)获得成功应用的重要原因之一在于核范数和L₁模的迫近算子均具有显式解。L₁模的迫近算子为众所周知的软阈值算子，设$Z\text{=pro}{{\text{x}}_{\beta {{\left\| \cdot \right\|}_{1}}}}\left(X \right)$，则其第(i, j)个元素可由下式显式给出：

$ {Z_{ij}} = \max \left\{ {\left| {{X_{ij}}} \right| - \beta ,0} \right\}{\rm{sign}}\left( {{X_{ij}}} \right)。$

(4)

注意到${{\left\| X \right\|}_{*}}={{\left\| \sigma \left(X \right) \right\|}_{1}}$，其中σ(X)表示X按降序排列的奇异值向量，文献^[11]中给出了核范数的迫近算子的显式解。设X=UΣV^T是X的奇异值分解。其中Σ=diag(σ(X))，则

$ {\rm{pro}}{{\rm{x}}_{\beta {{\left\| \cdot \right\|}_1}}}\left( X \right) = U\mathit{\tilde \Sigma }{V^{\rm{T}}}。$

(5)

其中

$ {{\mathit{\tilde \Sigma }}_{ii}} = \max \left\{ {\left| {{\mathit{\Sigma }_{ii}}} \right| - \beta ,0} \right\}。$

(6)

L₁模型虽然数值求解容易，但是对大于阈值β的元素将会引入偏差。文献[10]针对压缩感知问题改进了L₁模的上述缺点，提出了一种新的非凸稀疏罚函数：

$ {h_{\rho ,\tau }}\left( x \right) = \left\{ \begin{array}{l} \left( {2\rho - \tau } \right)\left| x \right| + {\left( {\rho - \tau } \right)^2},\left| x \right| \ge 2\left( {\tau - \rho } \right)\\ - \frac{1}{4}{x^2} + \rho \left| x \right|,\;\;\;\;\left| x \right| < 2\left( {\tau - \rho } \right) \end{array} \right.。$

(7)

其中ρ和τ满足$0 < \frac{\tau }{2}\le \rho \le \tau $。

从图 1可以看出h_ρ,τ(x)是L₀模的一种非凸近似函数并且其逼近效果优于L₁模。下面的定理1指出h_ρ,τ(x)的迫近算子具有显式解。

定理1  罚函数h_ρ,τ(x)的迫近算子可按照如下情形显式计算：

若β < 2，则prox_βhρ,τ(x)表达式如下

$ \left\{ \begin{array}{l} a - \beta \left( {2\rho - \tau } \right),a \ge \left( {2 - \beta } \right)\tau + \left( {2 - \beta } \right)\rho \\ \frac{2}{{2 - \beta }}\left( {a - \beta \rho } \right),\beta \rho < a < \left( {2 - \beta } \right)\tau + 2\left( {\beta - 1} \right)\rho \\ 0,\;\;\;\;\; - \beta \rho \le a \le \beta \rho \\ \frac{2}{{2 - \beta }}\left( {a + \beta \rho } \right), - \left( {2 - \beta } \right)\tau - 2\left( {\beta - 1} \right)\rho < a < - \beta \rho \\ a + \beta \left( {2\rho - \tau } \right),a \le - \left( {2 - \beta } \right)\tau - 2\left( {\beta - 1} \right)\rho \end{array} \right.。$

(8)

若β>2，则prox_βhρ,τ(x)表达式如下：

$ \left\{ \begin{array}{l} a - \beta \left( {2\rho - \tau } \right),a \ge \beta \rho \\ 2\left( {\tau - \rho } \right),\;\;\;\;\frac{{3\beta - 2}}{2}\rho + \frac{{2 - \beta }}{2}\tau < a < \beta \rho \\ 0,\; - \frac{{3\beta - 2}}{2}\rho - \frac{{2 - \beta }}{2}\tau \le a \le \frac{{3\beta - 2}}{2}\rho + \frac{{2 - \beta }}{2}\tau \\ 2\left( {\rho - \tau } \right), - \beta \rho < a < - \frac{{3\beta - 2}}{2}\rho - \frac{{2 - \beta }}{2}\tau \\ a + \beta \left( {2\rho - \tau } \right),a \le - \beta \rho \end{array} \right.。$

(9)

若β=2，则prox_βhρ,τ(x)表达式如下：

$ \left\{ \begin{array}{l} a - 2\left( {2\rho - \tau } \right),\;\;\;a \le \beta \rho \\ 0,\;\;\;\;\; - \beta \rho < a < \beta \rho \\ a + 2\left( {2\rho - \tau } \right),\;\;\;a \le - \beta \rho \end{array} \right.。$

(10)

证明 考虑到对称性，只需讨论a≥0的情形。由于h_ρ,τ(x)是分段函数，对应的迫近算子也是分段函数，分别用f(x)、g(x)代替。

令$f\left(x \right)=\frac{1}{2}{{\left(x-a \right)}^{2}}+\beta \left(2\rho-\tau \right)x+\beta {{\left(\rho-\tau \right)}^{2}}$，

$ \begin{array}{l} x \ge 2\left( {\tau - \rho } \right);g\left( x \right) = \frac{1}{2}{\left( {x - a} \right)^2} - \frac{\beta }{4}{x^2} + \beta \rho x,\\ 0 \le x < 2\left( {\tau - \rho } \right)。\end{array} $

对于f(x)，如a≥(2-β)τ+2(β-1)ρ，

$ \begin{array}{l} {f_{\min }} = f\left( {a - \beta \left( {2\rho - \tau } \right)} \right) = \\ - \frac{1}{2}\left( {2\rho - \tau } \right){\beta ^2} + \left[ {{{\left( {\rho - \tau } \right)}^2} + a\left( {2\rho - \tau } \right)} \right]\beta \end{array} $

如a < (2-β)τ+2(β-1)ρ，

$ \begin{array}{l} {f_{\min }} = f\left( {2\left( {\tau - \rho } \right)} \right) = \\ \left( {\tau - \rho } \right)\left[ {\left( {3\rho - \tau } \right)\beta + 2\left( {\tau - \rho - a} \right)} \right] + \frac{1}{2}{a^2}。\end{array} $

对于g(x)，我们分三种情况讨论g(x)在不同区间的最小值。

(1) β < 2：如βρ < a < (2-β)τ+2(β-1)ρ，

$ {g_{\min }} = g\left( {\frac{2}{{2 - \beta }}\left( {a - \beta \rho } \right)} \right) = - \frac{{{{\left( {a - \beta \rho } \right)}^2}}}{{2 - \beta }} + \frac{1}{2}{a^2} $

如a≥(2-β)τ+2(β-1)ρ，

$ \begin{array}{l} {g_{\min }} = g\left( {2\left( {\tau - \rho } \right)} \right) = \\ \left( {3\rho - \tau } \right)\left( {\tau - \rho } \right)\beta + 2\left( {\tau - \rho - a} \right)\left( {\tau - \rho } \right) + \frac{1}{2}{a^2} \end{array} $

如a≤βρ，${{g}_{\text{min}}}=g\left(0 \right)=\frac{1}{2}{{a}^{2}}$。

(2) β>2，如a≥βρ，

$ \begin{array}{l} {g_{\min }} = g\left( {2\left( {\tau - \rho } \right)} \right) = \\ \left( {3\rho - \tau } \right)\left( {\tau - \rho } \right)\beta + 2\left( {\tau - \rho - a} \right)\left( {\tau - \rho } \right) + \frac{1}{2}{a^2} \end{array} $

如$a < \left(2-\beta \right)\tau +2\left(\beta-1 \right)\rho $，${{g}_{\text{min}}}=g\left(0 \right)=\frac{1}{2}{{a}^{2}}$

如$\frac{3\beta-2}{2}\rho +\frac{2-\beta }{2}\tau < a < \beta \rho $，

$ \begin{array}{l} {g_{\min }} = g\left( {2\left( {\tau - \rho } \right)} \right) = \\ \left( {3\rho - \tau } \right)\left( {\tau - \rho } \right)\beta + 2\left( {\tau - \rho - a} \right)\left( {\tau - \rho } \right) + \frac{1}{2}{a^2} \end{array} $

如$\left(2-\beta \right)\tau +2\left(\beta-1 \right)\rho ~ < a < \frac{3\beta-2}{2}\rho +\frac{2-\beta }{2}\tau $，

$ {g_{\min }} = g\left( 0 \right) = \frac{1}{2}{a^2} $

(3) 当β=2时，

$ g\left( x \right) = \left( {2\rho - a} \right)x + \frac{1}{2}{a^2},0 \le x < 2\left( {\tau - \rho } \right) $

如a≤2ρ，${{g}_{\text{min}}}=g\left(0 \right)=\frac{1}{2}{{a}^{2}}$

如a>2ρ，${{g}_{\text{min}}}=g\left(2\left(\tau-\rho \right) \right)=$

$ 2\left( {\tau - \rho } \right)\left( {2\rho - a} \right) + \frac{1}{2}{a^2} $

综合上述分析，比较f_min和g_min的值，可得结果。

图 2画出L₁模、L₀模及β=1，ρ=0.75τ时的h_ρ,τ(x)对应的迫近算子的图像，可以看出h_ρ,τ(x)可以修正L₁模导致的对于大于阈值的变量的偏差。

基于h_ρ,τ(x)的上述性质，我们将其推广到矩阵或向量型数据，定义${{G}_{\rho, \tau }}\left(X \right)=\underset{x}{\mathop{\Sigma }}\, {{h}_{\rho, \tau }}\left(x \right)$，其中x为X的所有分量，提出如下的稀疏低秩矩阵分解的非凸模型：

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{L,S} {G_{\rho ,\tau }}\left( {\sigma \left( L \right)} \right) + \lambda {G_{\rho ,\tau }}\left( S \right)}\\ {{\rm{s}}.\;{\rm{t}}.\;\;\;L + S = D。} \end{array} $

(11)

1.2 基于ADMM的数值算法

ADMM是一类用于求解线性等式约束凸优化问题的有效算法，ADMM的思想可推广到求解非凸优化问题，且其收敛性已经有了部分结果^[12]。

考虑新模型(11)的增广Lagrange乘子：

$ \begin{array}{*{20}{c}} {{L_\beta }\left( {L,S,\mathit{\Lambda }} \right) = {G_{\rho ,\tau }}\left( {\sigma \left( L \right)} \right) + \lambda {G_{\rho ,\tau }}\left( S \right) - }\\ {\left[ {\mathit{\Lambda },L + S - D} \right] + \frac{\beta }{2}\left\| {L + S - D} \right\|_F^2。} \end{array} $

(12)

其中Λ∈R^m×n为线性约束乘子，β>0为惩罚因子。求解模型(11)的ADMM算法的迭代格式为：

$ \left\{ \begin{array}{l} {S^{k + 1}} = \mathop {\arg \min }\limits_S {L_\beta }\left( {S,{L^k},{\mathit{\Lambda }^k}} \right)\\ {L^{k + 1}} = \mathop {\arg \min }\limits_L {L_\beta }\left( {{S^{k + 1}},L,{\mathit{\Lambda }^k}} \right)\\ {\mathit{\Lambda }^{k + 1}} = {\mathit{\Lambda }^k} - \beta \left( {{L^{k + 1}} + {S^{k + 1}} - D} \right) \end{array} \right.。$

(13)

通过对更新L和S所需求解的子问题的目标函数进行分解与重组，并且借用迫近算子的定义，可得到式(13)的具体迭代形式。其中S^k+1的第(i, j)个元素可由下式计算：

$ S_{ij}^{k + 1} = {\rm{pro}}{{\rm{x}}_{\frac{\lambda }{\beta }{h_{\rho ,\tau }}}}\left( {{D_{ij}} - L_{ij}^k + \frac{1}{\beta }\mathit{\Lambda }_{ij}^k} \right)。$

(14)

设U^k+1Σ^k+1(V^k+1)^T是$D-{{S}^{k+1}}+\frac{1}{\beta }{{\mathit{\Lambda} }^{k}}$的奇异值分解，则L^k+1可由下式计算：

$ {L^{k + 1}} = {U^{k + 1}}{{\mathit{\tilde \Sigma }}^{k + 1}}{\left( {{V^{k + 1}}} \right)^T}。$

(15)

其中${{{\tilde{\Sigma }}}_{ii}}^{k+1}=\text{pro}{{\text{x}}_{\frac{1}{\beta }{{h}_{\rho, \tau }}}}({{\Sigma }_{ii}}^{k+1})$

Λ^k+1可由下式计算：

$ {\Lambda ^{k + 1}} = {\Lambda ^k} - \beta \left( {{L^{k + 1}} + {S^{k + 1}} - D} \right)。$

(16)

算法的运算效率取决于式(13)中的前两个子问题的求解。由于h_ρ,τ(x)的迫近算子具有显式解，式(14)、(15)中对S^k+1和L^k+1的更新均可以按照定理1中的结果显式给出。

2 数值实验

为了验证模型和算法的有效性，本节进行数值模拟实验。采用relerr指标来评估恢复结果

$ {\rm{relerr}} = \frac{{{{\left\| {\left( {\hat L,\hat S} \right) - \left( {\tilde L,\tilde S} \right)} \right\|}_F}}}{{{{\left\| {\left( {\tilde L,\tilde S} \right)} \right\|}_F}}}。$

其中$\left(\tilde{L}, \tilde{S} \right)$为原始矩阵，$\left(\hat{L}, \hat{S} \right)$为运行算法所恢复的矩阵。采用r和s表示模拟矩阵的秩和稀疏度，取m=n=100，本文基于如下的MATLAB命令生成：

$ D = \tilde L + \tilde S: $

$ \tilde L = {\rm{randn}}\left( {m,r} \right) \times {\rm{randn}}\left( {r,n} \right); $

$ S = {\rm{zeros}}\left( {m,n} \right);p = {\rm{randperm}}\left( {m \times n} \right); $

$ S' = {\rm{round}}\left( {s \times m \times n} \right); $

$ \tilde S\left( {p\left( {1;S'} \right)} \right) = \max \left( {abs\left( {\tilde L\left( : \right)} \right)} \right) \times {\rm{sign}}\left( {{\rm{randn}}\left( {S',1} \right)} \right); $

正则化参数λ对模型恢复结果具有显著影响，但是目前尚缺乏对该参数的有效的自适应选取策略，本文通过数值实验手动选取能够产生最优分解的正则化参数。参考文献[6]中的方法取$\beta =0.25\text{mn}/{{\left\| D \right\|}_{1}}$，当相邻两次迭代误差满足如下停机准则时，算法终止

$ \frac{{{{\left\| {\left( {{L^{k + 1}},{S^{k + 1}}} \right) - \left( {{L^k},{S^k}} \right)} \right\|}_F}}}{{{{\left\| {\left( {{L^k},{S^k}} \right)} \right\|}_F}}} \le {10^{ - 6}}。$

2.1 参数ρ、τ的选取

下面讨论参数ρ、τ在不同取值下对模型恢复能力的影响。由于ρ∈[τ/2, τ]，令σ=ρ/τ，图 3描述了针对几种模拟实验，σ在区间[1/2, 1]内的不同取值下恢复误差relerr的变化情况。可以看出，当σ>0.6时，其值对模型的恢复结果没有显著影响。本文下面实验中选取σ=0.75。

2.2 模型分解能力的比较

下面比较本文模型(11)和广泛使用的L₁模型(2)的恢复结果，其中L₁模型采用文献[6]中的算法求解。为充分对比不同模型的恢复能力，我们模拟它们的恢复成功率。设置r在5~50之间变化，s在5%~40%之间变化，对于每一对(r, s)，重复数值实验50次并且计算恢复误差，当relerr < 10^-3时，认为恢复矩阵L和S成功。令p表示恢复成功的次数，则恢复成功率定义为p/50。图 4和5分别显示了L₁模型(2)和本文模型(11)的恢复成功率。从图 4和5可以看出，当r较小或者是s较小时，两模型均可以成功地恢复矩阵的低秩和稀疏成分。当稀疏度较小时，本文模型对低秩矩阵的恢复能力更强。例如，当s=0.1时，本文模型可成功恢复r≤30时的低秩矩阵，而L₁模型只能成功恢复r < 25时的低秩矩阵。然而，当稀疏度s较大时，本文模型的恢复成功率略差于L₁模型。

3 结语

本文针对稀疏低秩矩阵分解问题提出了一种新的非凸近似模型，并且基于ADMM算法发展了模型的数值解法。数值实验结果表明，本文模型相较于广泛使用的L₁凸近似模型具有一定的优势。

考虑到模型中的目标函数非凸，本文数值算法的理论收敛性等结果尚未完全建立，模型中相关参数的选择亦缺乏相应的自适应策略。在以后的研究中，我们将研究模型的理论恢复能力，发展可靠、高效的数值算法，并且建立参数的自适应选取策略，以使得模型可以方便地解决稀疏低秩矩阵分解的实际问题。