有效地解决低阶矩阵完全问题（对矩阵的ky-norm做惩罚）
Solving Low-Rank Matrix Completion Problems Efficiently
Donald Goldfarb Shiqian Ma（马士谦，香港中文大学） Zaiwen Wen（文再文，北京大学）
摘要：我们提出(present)几个求解低秩矩阵完全问题的一阶算法，通过最小化核范数(nuclear norm)代替矩阵的秩得到的它的最紧凸松弛(tightest convex relaxation)。我们的第一个算法是一个不动点延拓算法(fixed point continuation algorithm)，它整合了一个近似奇异值分解过程(FPCA, approximate singular value decomposition procedure)。FPCA可以有效地解决大型矩阵完全问题，并达到高水平的可恢复性(recoverability)。例如，FPCA可以在大约3分钟内恢复秩为50的$1000\times 1000$矩阵，相对误差为$10^{-5}$，仅仅需要对$20\%$的元素进行抽样。我们知道没有其他的方法能够达到这么良好的可恢复性。我们的第二个算法是求解核范数矩阵完全问题的半定规划重写(semidefinite programming reformulation)的逐行法(row by row method)。这种方法能够高效地产生相当大的核范数矩阵完全问题的高精度解(highly accurate solutions)。最后，我们介绍一种基于增广拉格朗日框架(augmented Lagrangian framework)的交替方向方法。
1 介绍
本文中，我们对求解矩阵完全问题的方法感兴趣。
\begin{align*}&\min &&\rank (X)\\
&s.t. &&X_{ij}=M_{ij},\quad \forall (i,j)\in \Omega,
\end{align*}
其中$X$和$M$都是$p\times q$矩阵和$\Omega$是指标对$(i,j)$的子集.此问题是仿射约束矩阵秩最小化问题(affinely constrained matrix rank minimization problem)
\begin{align*}&\min &&\rank (X)\\
&s.t. &&\mathcal{A}(X)=b,
\end{align*}
的特殊情形，其中$X\in \mathbb{R}^{p\times q}$为决策变量，而线性映射$\mathcal{A}: \mathbb{R}^{p\times q}\to \mathbb{R}^m$且向量$b\in \mathbb{R}^m$是给定的.
所谓的协同过滤问题(collaborative filtering problem) [17]可以被视为(be cast as)一个矩阵完全问题.假设用户在一个在线调查提供一些电影的评分.这产生一个以用户为行和电影为列的矩阵$M$，它的$(i,j)$项$M_{ ij}$是第$i$位用户对第$j$部电影给出的评分.因为大多数用户只对一小部分的电影评分，我们通常只知道矩阵所包含项的一个小的子集$\{M_{ij}|(i,j)\in \Omega\}$.根据某个用户的已知评分，我们想要预测用户对那些他没评过分的电影的评分；也就是，我们要填写矩阵的那些空白项.通常认为，只有少量因素影响一个人对电影的品味或偏好.因此评分矩阵$M$可能从数值上看是低秩的，在此意义上，相对较少的大奇异值占(account for)所有的奇异值的总和的很大一部分。找到这样一个低秩矩阵$M$对应于解决矩阵完全问题(I.1).
矩阵的秩是其正奇异值的个数.矩阵秩最小化问题(I.2)一般是NP难的，这是由函数$\rank (\cdot)$的组合性质导致的。为了得到一个便于计算的近似(I.2),我们可以用它的凸包来代替$\rank (X)$,也就是核范数$\|X\|_\ast$,它被定义为矩阵$X$所有奇异值的总和 [8].由此产生了以下核范数最小化问题，这是(I.2)的最紧凸松弛：
\begin{align*}&\min &&\|X\|_\ast\\
&s.t. &&\mathcal{A}(X)=b.
\end{align*}
如果$b$被噪声污染(contaminated by)了，约束(constraint) $\mathcal{A}(X)=b$应该被放宽，也就产生了另一个问题
\begin{align*}&\min &&\|X\|_\ast\\
&s.t. &&\|\mathcal{A}(X)-b\|_2\leq \theta.
\end{align*}
或者是拉格朗日函数的说法(version)
\[\min \mu\|X\|_\ast+\frac12\|\mathcal{A}(X)-b\|_2^2,\]
其中$\theta$和$\mu$是参数.
一个重要的问题是：核范数最小化问题(I.3)的最优解啥时候能给出矩阵秩最小化问题(I.2)的最优解.针对这一问题，Recht等人[15]证明了：若$A$ ($A$为线性算子$\mathcal{A}$的矩阵说法，也就是$\mathcal{A}(X)=A\mathrm{vec} (X)$)的各项是适当随机的(suitably random), 例如，独立同分布的正态分布，那么只要$m\geq Cr (p+q)$,大部分秩为$r$的$p\times q$矩阵均可被复原为求解核范数最小化问题(I.3),其中$C$为正常数.
对于矩阵完全问题(I.1),相应的核范数优化问题为
\begin{align*}&\min &&\|X\|_\ast\\
&s.t. &&X_{ij}=M_{ij},\forall (i,j)\in \Omega.
\end{align*}
Candès和Tao在[6]证明了在某些无序条件下(under certain incoherence conditions),低秩矩阵完全可以以高概率复原为核范数最小化问题(I.5),如果样本数$m$的阶为$O(nr\mathrm{polylog}(n))$,其中$n=\max (p,q)$.
在本文的其余部分，我们提出一些高效的方法来求解(I.1), (I.2), (I.3), (I.4)和 (I.5).
II 不动点迭代算法
我们求解(I.4)的不动点迭代算法是如下简单的两行算法
\[
\begin{cases}
Y^k=X^k-\tau g(X^k)\\
X^{k+1}=S_{\tau \mu}(Y^k),
\end{cases}
\]
其中$g(X^k)=\mathcal{A}^\ast (\mathcal{A}(X^k)-b)$, $\mathcal{A}^\ast$为$\mathcal{A}$的伴随算子，且$S_\nu (\cdot)$为矩阵收缩算子(matrix shrinkage operator),定义为
\[S_\nu (Y)=U\mathrm{Diag}(\bar\sigma)V^T,\quad \bar\sigma=\max\{\sigma-\nu,0\},\]
其中$ U\mathrm{Diag}(\bar\sigma)V^T$是矩阵$Y\in \mathbb{R}^{p\times q}$的奇异值分解(SVD).
我们的算法(II.1)灵感来源于[12]提出的不动点迭代算法。因为(I.4)中的目标函数(objective function)是凸的，$X^\ast$是(I.4)的最优解(optimal solution)当且仅当
\[0\in \mu\partial\|X^\ast\|_\ast+g(X^\ast).\]
注意到如果$X$的SVD为$X=U\Sigma V^T$,则(见[3])
$$
\partial \lVert X \rVert _{\ast}=\left\{ UV^T+W\left| U^TW=0,WV=0,\lVert W \rVert _2\leq 1 \right. \right\},
$$
其中范数$\|W\|_2$定义为$W$的最大奇异值.
基于最优条件(II.3),我们可以通过采用算子分解技巧(an operator splitting technique)开发出一种求解(I.4)的不动点迭代方法.注意到(II.3)等价于对任意$\tau>0$,均有
\[0\in \tau\mu\partial \|X^\ast\|_\ast+X^\ast-( X^\ast-\tau g (X^\ast)).\]
若令\[Y^\ast=X^\ast-\tau g(X^\ast),\]则(II.4)变为
\[0\in \tau\mu\partial \|X^\ast\|_\ast+X^\ast-Y^\ast,\]
即$X^\ast$是
\[\min \tau\mu \|X^\ast\|_\ast +\frac12 \left\|X -Y^\ast\right\|_\mathrm{F}^2. \]
我们知道应用于$Y^\ast$的矩阵收缩算子，也就是$S_{\tau\mu}(Y^\ast)$,给出了(II.7)的最优解.（见[13]）
基于以上结论我们得到本节一开始提到的求解问题(I.4)的不动点迭代算法(II.1).而且，我们有以下定理：
定理2.1. $X^\ast$是问题(I.4)的最优解当且仅当$X^\ast=S_{\tau\mu}(h(X^\ast))$,其中$h(\cdot)=I(\cdot)-\tau g(\cdot)$.
- 收敛性结果
我们现在分析不动点迭代算法的收敛性（证明见[13]）.
引理2.2. 收缩算子$S_\nu$是非扩张的(nonexpansive),也就是对任何$Y_1$和$Y_2\in \mathbb{R}^{p\times q},\nu>0$,均有
\[\left\|S_\nu (Y_1)-S_\nu (Y_2)\right\|_{\mathrm{F}}\leq \| Y_1- Y_2\right\|_{\mathrm{F}}.\]
而且，
$$
\lVert Y_1-Y_2 \rVert _F=\lVert S_{\nu}\left( Y_1 \right) -S_{\nu}\left( Y_2 \right) \rVert _F\Leftrightarrow Y_1-Y_2=S_{\nu}\left( Y_1 \right) -S_{\nu}\left( Y_2 \right) .
$$
我们现在说明不动点迭代(II.1)收敛到问题(I.4)的最优解.
定理2.3. 由不动点迭代算法(II.1)生成的数列$\left\{ X^k \right\}$收敛到某些$X^\ast \in\mathcal{X}^\ast$,其中$\mathcal{X}^\ast$是问题(I.4)最优解集,而$\tau\in (0,2/\lambda_\max (\mathcal{A}^\ast\mathcal{A}))$.
- 延拓性(Continuation)
受Hale [12]工作的启发，我们首先描述一个延拓（同伦）(continuation (homotopy))技巧，以加快不动点迭代(II.1)的收敛.我们求解(I.4)的不动点延拓(FPC, fixed point continuation)迭代算法如下所示：
- 基于FPC算法的一种近似SVD：FPCA
计算奇异值分解是算法1的主要计算代价.因此，为了取代在每次迭代中计算矩阵$Y$的完整SVD,我们执行了算法FPC的一个变体，这样的话我们只计算秩为$k_s$的$Y$的近似（见下文），其中$k_s$是一个在算法过程中调整的参数.我们称其为变体FPCA.
这种方法大大降低了算法所需的计算量.具体地说，我们通过快速蒙特卡罗算法(fast Monte Carlo algorithm)计算一个近似的SVD：Drineas等[7]所开发的线性时间SVD算法.对于一个给定的矩阵$A\in \mathbb{R}^{p\times q}$,和参数$c_s,k_s\in\mathbb{Z}^+$,其中$1\leq k_s\leq c_s\leq q$且$\{p_i\}_{i=1}^q,p_i\geq 0,\sum_{i=1}^q p_i=1$,该算法在线性$O (p+q)$时间内返回一个矩阵$A$的最大$k_s$奇异值和相应的左奇异向量的近似.
下面概述了线性时间SVD算法.