奇异值阈值算法

本文提出了一种新的算法来逼近所有服从凸约束的矩阵中具有最小核范数的矩阵。这个问题可以理解为一个秩最小化问题的凸松弛，并且在许多重要的应用中出现，比如从其条目(著名的Netflix问题)的一个小子集恢复一个大矩阵。像内点法这样的现成算法并不直接适用于这类含有上百万个未知项的大型问题。

现实生活中的大规模数据常常会有部分数据缺失、数据误差、损坏等问题，这将进一步加大数据处理和分析难度。这种现在在实际生活中很常见，例如在人脸识别中,训练集中的或是待识别的人脸图像往往含有阴影、高光、遮挡、变形等;在运动恢复结构(Structure from motion,SFM)问题中,进行特征提取和特征匹配时往往存在大的匹配误差. 这些因素的存在使得很多传统的分析和处理方式失效, 需要新的理论和实用的算法为相关的应用提供理论支撑和有力的求解工具.正确并高效地从不完整的、带有损毁的数据中 恢复和利用它们, 对现代大规模数据的分析与处理至关重要。

假定原始数据矩阵是低秩的，但是矩阵中含有很多未知的元素。从一个不完整的矩阵中恢复出一个完整的低秩阵，便是低秩矩阵填充问题。例如，著名的Netflix问题便是一个典型的低秩矩阵填充问题。Netflix是美国的一家影片租赁公司，其推荐系（recommendation system）要从用户仅有的对少数的电影打分中为用户推荐影片，如果这种推荐越符合用户的喜好，也就越能提高该公司租赁电影的业务，为此，该公司设立了百万美元的奖金用于悬赏能够最好地提高该公司推荐系统质量的解决方法。这个问题可以用矩阵填充来进行建假设矩阵的每一行分别代表同一用户对不同电影的打分, 每一列代表不同用户对同一电影的打分，用户数量巨大, 电影数目巨大，因此这个矩阵的维度十分大。由于用户所打分的电影有限，这个矩阵中只有很小一部分的元素值已知，而且可能含有噪声或误差，那么Netflix问题就是如何从这个不完整的矩阵中推测其中未知元素的值。矩阵填充的越准确，为用户推荐的电影也就越符合用户的喜好。

低秩矩阵补全

矩阵补全是推荐系统、计算机视觉、图像处理等领域经常遇到的问题，具有很高的研究价值。矩阵补全问题旨在通过真实未知矩阵 M的部分观测 $\mathbf{M}_{i,j},(i,j)\in\Omega $​来估计 M 中未观测到的其他元素。如果不加其他约束，这样的问题是完全不可解的。但是，如果数据矩阵具有一些特殊的性质，这些特殊的性质将使得矩阵补全问题成为可能。低秩性就是这样的性质。如果事先知道矩阵 M是低秩的，那么矩阵补全问题可以形式为为如下优化问题：
KaTeX parse error: Can't use function '$' in math mode at position 2: $̲\begin{aligned}…
其中 其中 Ω 为观测样本下标的集合，X为优化变量，M为真实的未知矩阵。定义投影算子 $P_{\Omega }$:
$${\left[P_{\Omega }(\mathbf{X})\right]}_{i,j}=\{\begin{array}{ll}{\mathbf{X}}_{i,j}& (i,j)\in \Omega \\ 0& \mathrm{otherwise}\end{array}(2)$$
从而优化问题（1）可以简洁的描述为：
$$\begin{array}{rl}\underset{\mathbf{X}}{\min }& rank(\mathbf{X})\\ s.t.& P_{\Omega }(\mathbf{X})=P_{\Omega }(\mathbf{M})\end{array}(3)$$
显然，上述问题是一个NP-hard问题，且其复杂度随矩阵维度的增城呈平方倍指数关系增长。为了求解问题(1)，有学者对其进行凸松弛，转化为一个凸优化问题：
$$\begin{array}{rl}& \underset{\mathbf{X}}{\min }{\Vert \mathbf{X}\Vert }_{\ast }\\ & s.t.{\mathbf{X}}_{i,j}={\mathbf{A}}_{i,j},(i,j)\in \Omega \end{array}(4)$$
问题（4）的解在满足强非相干性条件下，很大概率等价于问题(1)的解。在某种意义上，这是NP-难秩最小化问题的紧凸松弛。因为核范数球是谱范数为1的rank one矩阵集合的凸包。具体定理及证明参见文献[1]。
对于(4)中的优化问题，当采样点数满足 $m\ge C{n}^{6/5}r\log n$ 时，M 能以很大概率求得精确解，其中 r为矩阵 M的秩，C 为一个正常数。

求解(4)式可以使用一些凸松弛方法，比如半定规划法，然而半定规划法通常使用内点法来求解，其求解上述问题的算法复杂度为 $O(p(m+n{)}^3+p^2(m+n{)}^2+p^3)$。可以通过matlab软件包CVX等来求解

奇异值阈值算法

奇异值收缩算子

考虑到秩为r的矩阵$X\in {\mathbb{R}}^{n_1\times n_2}$的奇异值分解：
X = U Σ V ∗ , Σ = d i a g ( { σ i } 1 ≤ i ≤ r ) X=U\Sigma V^*,\quad\Sigma=\mathrm{diag}(\{\sigma_i\}_{1\leq i\leq r}) X=UΣV∗,Σ=diag({σi​}1≤i≤r​)
其中U和V分别是具有正交列的n1 × r和n2 × r矩阵，奇异值σ i为正(除非另有说明,我们将始终假设矩阵的奇异值分解为上式的简化形式) .对于每一个τ≥0，我们引入如下定义的软阈值算子${\mathcal{D}}_{\tau }$​：
D τ ( X ) : = U D τ ( Σ ) V ∗ , D τ ( Σ ) = d i a g ( { σ i − τ ) + } ) \mathcal{D}_{\tau}(X):=U\mathcal{D}_{\tau}(\Sigma)V^{*},\quad\mathcal{D}_{\tau}(\Sigma)=\mathrm{diag}(\{\sigma_{i}-\tau)_{+}\}) Dτ​(X):=UDτ​(Σ)V∗,Dτ​(Σ)=diag({σi​−τ)+​})
其中t +是t的正部分，即$t_{+}=\max (0,t)$。换句话说，该算子简单地将软阈值规则应用于X的奇异值，有效地将这些奇异值收缩到零。这也是我们将这种变换称为奇异值收缩算子的原因。尽管奇异值分解可能并不唯一，但不难看出奇异值收缩算子的定义是很好的，我们不再赘述。从某种意义上说，这个收缩算子是对标量和向量的软阈值规则的直接扩展。特别地，注意到如果X的许多奇异值低于阈值τ，那么${\mathcal{D}}_{\tau }$的秩可能比X的秩要低得多，就像对向量施加软阈值规则会导致当输入的某些项低于阈值时输出更稀疏一样。

奇异值阈值算子是与核范数相关的邻近算子。关于邻近算子的详细信息。

**定理：**对于每个$\tau \ge 0$和$Y\in {\mathbb{R}}^{n_1\times n_2}$，奇异值收缩算子 D τ ( X ) : = U D τ ( Σ ) V ∗ , D τ ( Σ ) = d i a g ( { σ i − τ ) + } ) \mathcal{D}_{\tau}(X):=U\mathcal{D}_{\tau}(\Sigma)V^{*},\quad\mathcal{D}_{\tau}(\Sigma)=\mathrm{diag}(\{\sigma_{i}-\tau)_{+}\}) Dτ​(X):=UDτ​(Σ)V∗,Dτ​(Σ)=diag({σi​−τ)+​})服从
$${\mathcal{D}}_{\tau }(\mathbf{Y})=\arg \underset{\mathbf{X}}{\min }\left\{\frac{1}{2}\Vert \mathbf{X}-\mathbf{Y}{\Vert }_F^2+\tau \Vert \mathbf{X}{\Vert }_{\ast }\right\}.$$
证明：

由于函数$h_0(X):=\tau \Vert X{\Vert }_{\ast }+\frac{1}{2}\Vert X-Y{\Vert }_F^2$是严格凸的，容易看出存在唯一的极小点，因此需要证明它等于${\mathcal{D}}_{\tau }(Y)$。为此，回想一下凸函数f的次梯度定义：$f:{\mathbb{R}}^{n_1\times n_2}\to \mathbb{R}.$。称Z是f在$X_0$处的次梯度，记$Z\in \partial f(X_0)$，如果：
$$f(X)\ge f(X_0)+⟨Z,X-X_0⟩$$
对所有的X，现在$\hat{X}$极小化$h_0$当且仅当0是泛函h0在点$\hat{X}$的次梯度，即.

$$0\in \hat{X}-Y+\tau \partial \Vert \hat{X}{\Vert }_{\ast }$$
式中$\partial \Vert \hat{X}{\Vert }_{\ast }$为核范数的次梯度集合。设$X\in {\mathbb{R}}^{n_1\times n_2}$为任意矩阵，$U\Sigma V^{\ast }$为其SVD。众所周知，
$$\partial \Vert X{\Vert }_{\ast }=\left\{\begin{array}{ccccc}U{V}^{\ast }+W:& W\in {\mathbb{R}}^{n_1\times n_2},& U^{\ast }W=0,& WV=0,& \Vert W{\Vert }_2\le 1\end{array}\right\}.$$

记$\hat{X}:=D_{\tau }(Y)$ .为了说明$\hat{X}$服从$0\in \hat{X}-Y+\tau \partial \Vert \hat{X}{\Vert }_{\ast }$​，对Y进行SVD分解
$$Y=U_0{\Sigma }_0{V}_0^{\ast }+U_1{\Sigma }_1{V}_1^{\ast },$$
式中：$U_0,V_0(resp.U_1,V_1)$为奇异值大于 $\tau$(和小于或等于$\tau$)的奇异向量。有了这些记号，我们就有了

$$\hat{X}=U_0({\Sigma }_0-\tau I)V_0^{\ast }$$
因此
$$Y-\hat{X}=\tau (U_0{V}_0^{\ast }+W),W={\tau }^{-1}{U}_1{\Sigma }_1{V}_1^{\ast }.$$
由定义可知，$U_0^{\ast }W=0,W{V}_0=0$，又因为${\Sigma }_1$的对角元素有以τ为界的量值，所以$\Vert W{\Vert }_2\le 1$。从而$Y-\hat{X}\in \tau \partial \Vert \hat{X}{\Vert }_{\ast }$，证明结束.

奇异值阈值算法（SVT）

该算法可以类比为求解向量 ${\ell }_0$​范数最小化的软阈值迭代算法。 SVT算法先将最优化问题(4)正则化，即有：
$$\begin{array}{rl}\underset{\mathbf{X}}{\min }& \tau {\Vert \mathbf{X}\Vert }_{\ast }+\frac{1}{2}{\Vert \mathbf{X}\Vert }_F^2\\ s.t.& P_{\Omega }(\mathbf{X})=P_{\Omega }(\mathbf{M})\end{array}(5)$$
其中，$\tau >0$。当$\tau\to+\infty $​ 时，上述最优化问题的最优解收敛到(4)式的最优解。构造最优化问题(5)的拉格朗日函数：
$$L(\mathbf{X},\mathbf{Y})={\Vert \mathbf{X}\Vert }_{\ast }+\frac{1}{2}{\Vert \mathbf{X}\Vert }_F^2+⟨\mathbf{Y},P_{\Omega }(\mathbf{M}-\mathbf{X})⟩(6)$$
其中，拉格朗日乘子$\mathbf{Y}\in {\mathbb{R}}^{m\times n}$。如果 $({\mathbf{X}}^{\ast },{\mathbf{Y}}^{\ast })$​为优化问题(4)的原-对偶问题的最优解，则有：
$$\underset{\mathbf{Y}}{\sup }\underset{\mathbf{X}}{\inf }L(\mathbf{X},\mathbf{Y})=\underset{\mathbf{X}}{\inf }\underset{\mathbf{Y}}{\sup }L(\mathbf{X},\mathbf{Y})=L({\mathbf{X}}^{\ast },{\mathbf{Y}}^{\ast })(7)$$
SVT算法使用交替迭代方法求解优化问题(4)，其迭代格式可以简单表述如下:
$$\{\begin{array}{c}{\mathbf{X}}^k=D_{\tau }({\mathbf{Y}}^{k-1})\\ {\mathbf{Y}}^k={\mathbf{Y}}^{k-1}+{\delta }_k{P}_{\Omega }(\mathbf{M}-{\mathbf{X}}^k)\end{array}(8)$$
其中$D_{\tau }(\mathbf{W})$ 为奇异值阈值软阈值操作类似，不过这里的对象是矩阵。可以具体描述为
$${\mathbf{W}}^{\prime }=D_{\tau }(\mathbf{W})=\{\begin{array}{cc}(1)& [\mathbf{U},\mathbf{S},\mathbf{V}]=\mathrm{svd}(W)\\ (2)& \mathbf{S}=\mathrm{sgn}(\mathbf{S}).\ast \max (\mathrm{abs}(\mathbf{S})-\tau ,0)\\ (3)& {\mathbf{W}}^{\prime }=\mathbf{U}\ast \mathbf{S}\ast {\mathbf{V}}^T\end{array}(9)$$
其中(8)式中的第一步是这样来的，初始化${\mathbf{Y}}^0=0$，当 ${\mathbf{Y}}^{k-1}$​固定时：
$$\begin{array}{rl}{\mathbf{x}}^k& =\underset{\mathbf{X}}{\mathrm{argmin}}\tau {\Vert \mathbf{X}\Vert }_{\ast }+\frac{1}{2}{\Vert \mathbf{X}\Vert }_F^2-⟨{\mathbf{Y}}^{k-1},P_{\Omega }(\mathbf{X})⟩\\ & =\underset{\mathbf{X}}{\mathrm{argmin}}\tau \Vert \mathbf{X}{\Vert }_{\ast }+\frac{1}{2}\Vert \mathbf{X}{\Vert }_F^2-⟨\mathbf{X},P_{\Omega }({\mathbf{Y}}^{k-1})⟩(10)\\ & =D_{\tau }({\mathbf{Y}}^{k-1})\end{array}$$
(8)式中的第(2)式，当 ${\mathbf{x}}^k$ 给定时，用梯度下降来更新 Y。

ps:

邻近算子

定义邻近算子：
$${\mathrm{prox}}_h(x)=\underset{u}{\mathrm{argmin}}\left(h(u)+\frac{1}{2}\Vert u-x{\Vert }_2^2\right)$$
直观理解：求解一个距x不算太远的点u，并使函数值h(u)也相对较小

数学理解：对下列最小化问题
$$\underset{\mathbf{u}}{\min }\frac{1}{2}{\Vert \mathbf{u}-\mathbf{x}\Vert }_2^2+h(\mathbf{u})$$
存在唯⼀的解$\hat{u}$ ，并将唯⼀的最小点记为 ${\mathrm{prox}}_h(x)$。称 ${\mathrm{prox}}_h(x)$为函数 h(u) 在点 $:x\in {\mathbb{R}}^n$处的邻近算⼦。也称 ${\mathrm{prox}}_h(x)$是点 x关于函数 h(u) 的邻近点。

定理：（邻近算子是良定义的）

如果 h 为闭凸函数，则对任意 x，${\mathrm{prox}}_h(x)$​存在且唯一

证明

首先注意到$h(u)+\frac{1}{2}\Vert u-x{\Vert }_2^2$​ 是强凸函数，则

• 存在性：强凸函数的所有α-下水平集有界，故由Weierstrass定理 知最小值存在
• 唯一性：强凸函数最小值唯一

邻近算子与次梯度的关系

定理：

若h 是适当的闭凸函数, 则$ u =\mathrm{prox}_h(x)⇐⇒ x − u ∈ ∂h(u)$​

证明

若$u={\mathrm{prox}}_h(x)$, 则由最优性条件得0 ∈ ∂h(u) + (u − x)(存在某个 𝑔∈∂ℎ(𝑢) 使得0=𝑔+(𝑢−𝑥)，或等价地，−𝑔=𝑢−𝑥，或者𝑢=𝑥+𝑔,这种条件可以视为一种最优性条件，表明 𝑢 是某种平衡或均衡点。)

因此有$x-u\in \partial h(u)$.反之，若$x-u\in \partial h(u)$​，则由次梯度定义可得
$$h(\nu )⩾h(u)+(x-u{)}^{\mathrm{T}}(\nu -u),\forall \nu \in \mathrm{dom}h$$
两边同时加上$\frac{1}{2}\Vert \nu -x{\Vert }^2$，即有
$$\begin{array}{c}h(\nu )+\frac{1}{2}\Vert \nu -x{\Vert }^2⩾h(u)+(x-u{)}^{\mathrm{T}}(v-u)+\frac{1}{2}\Vert (v-u)-(x-u){\Vert }^2\\ ⩾h(u)+\frac{1}{2}\Vert u-x{\Vert }^2,\forall \nu \in \text{dom }h\end{array}$$
(因为$x-u\in \partial h(u)$所以由泰勒公式的一阶展开可得$h(\nu )\approx h(u)+(x-u{)}^{\mathrm{T}}(\nu -u)$ 而$\frac{1}{2}\Vert (v-u)-(x-u){\Vert }^2\ge \frac{1}{2}\Vert u-x{\Vert }^2$表示平方距离的增量，因此整体的左边表达式大于等于右边的表达式。）

根据定义可得$u={\mathrm{prox}}_h(x)$​

次梯度

次梯度方法(subgradient method)是传统的梯度下降方法的拓展，用来处理不可导的凸函数。它的优势是比传统方法处理问题范围大，劣势是算法收敛速度慢。但是，它对不可导函数有很好的处理方法。

定义

对于一个凸函数$f:{\mathbb{R}}^n\to \mathbb{R},$，在某些点上它可能不光滑，即不可微。在这些点上，我们不能简单地定义一个唯一的梯度，而是可以定义一个集合，这个集合被称为该点的次梯度集合

具体来说，对于凸函数 𝑓，在点$x\in {\mathbb{R}}^n$*，其次梯度集合**$\partial f(x)$定义为包含所有满足以下条件的向量$g\in {\mathbb{R}}^n$ 的集合**：
$$f(y)\ge f(x)+g^T(y-x),\forall y\in {\mathbb{R}}^n$$
则称 𝑔 是函数 𝑓 在点 𝑥 处的次梯度（Subgradient）。换句话说，𝑔是一个满足上述不等式关系的向量。这个条件反映了 𝑔是一个沿着某个方向使得函数值增加的方向，它可以在凸函数的不可微点提供一种广义的梯度概念。

计算

设函数 𝑓 在点 𝑥0 处不一定可导，以一维情况为例，按高等数学中我们对导数的定义可以求出:

f在$x_0$处的左导数：
$$a=\underset{x\to x_0^{-}}{\lim }\frac{f(x)-f(x_0)}{x-x_0}$$
f在$x_0$处的右导数：
$$b=\underset{x\to x_0^{+}}{\lim }\frac{f(x)-f(x_0)}{x-x_0}$$
凸函数 𝑓 的次微分等于闭区间 [𝑎 , 𝑏] 中任何一个取值都是次梯度

如果凸函数 𝑓 在点 𝑥 处是可导的，即 𝑎=𝑏 ，次微分中只有一个元素，此时次梯度就是梯度，即 𝑔 就等于 ∇𝑓(𝑥) ;

如果凸函数 𝑓 在点 𝑥 处是不可导的，即 𝑎≠𝑏 ，此时次梯度是次微分中的任意一个取值，它是不唯一的