几种矩阵完备算法的研究与实现——《矩阵分析》课程仿真作业报告*刘鹏飞电⼦系2016210858摘要矩阵完备是指从⼀⼩部分已知的矩阵元素中恢复出整个矩阵。
它在计算机视觉、推荐系统以及社交⽹络等⽅⾯具有⼴泛的应⽤。
矩阵恢复可以通过求解⼀个与核范数有关的凸优化问题来实现。
由此诞⽣了许多矩阵恢复的算法，⽐如FPC算法等。
FPC算法虽然实现简单，但其迭代速度较慢。
在此基础上，APG算法经过改进，能够提升迭代速度。
但最⼩化核范数并不是求解矩阵完备问题的唯⼀⽅法，其中OptSpace算法构造了⼀个在流形上的优化问题，相⽐于前两种算法能够以更⾼的精度恢复出原始矩阵。
本⽂主要总结了FPC、APG和OptSpace三种算法的步骤。
特别地，对于OptSpace算法，本⽂提出了求解其中两个⼦优化问题的具体算法。
最后，本⽂通过仿真实验和理论分析⽐较了三种算法的特点，并给出了OptSpace算法的精度⾼于APG算法的解释。
关键词：矩阵完备，核范数，FPC，APG，OptSpace1介绍1.1矩阵完备及其算法综述矩阵完备是指从⼀⼩部分已知的矩阵元素中恢复出整个矩阵。
它在计算机视觉、推荐系统以及社交⽹络等⽅⾯具有⼴泛的应⽤。
矩阵完备可以描述成这样⼀个问题：对于⼀个m×n的矩阵M，其秩为r，我们只有对M中的部分采样，记*报告中所涉及到的仿真代码可在https:///s/1jHRcY8m下载1这些采样位置组成的集合为Ω，那么是否有可能从已知的部分元素中恢复出整个矩阵M。
假如M为低秩矩阵，并且已知的元素⾜够多并且⾜够均匀地分布在整个矩阵中，那么我们可以通过解如下优化问题来恢复出原始矩阵[1]：min rank(W)s.t.W ij=M ij,(i,j)∈Ω(1-1)但是，问题(1-1)是⼀个NP难的⾮凸问题。
在⼀定条件下，问题(1-1)可以转化成⼀个最⼩化核范数的问题。
对于矩阵W m×n，W的核范数定义为其奇异值之和，即∥W∥∗=min(m,n)∑k=1σk(W)(1-2)其中，σk(W)表⽰W第k⼤的奇异值。
问题(1-1)可以转化成：min∥W∥∗s.t.W ij=M ij,(i,j)∈Ω(1-3)对于(1-3)中带等式约束的问题，进⼀步地，可以将它凸松弛成⼀个⽆约束的优化问题[2][3][4]：min 12∥A(W)−b∥22+µ∥W∥∗(1-4)其中，b是由矩阵中采样位置对应的元素组成的p×1维向量，p=|Ω|（|·|表⽰集合的势）；A:R m×n−→R p是⼀个线性映射，A(W)=(W ij)|(i,j)∈Ω；µ是⼀个可以调整的参数。
对于(1-4)中的⽆约束问题，⽂献[2][3]分别提出了Fixed Point Continuation （FPC）和Singular Value Thresholding（SVT）的算法。
本⽂认为，这两种算法虽然出发点不同，但其实质都是梯度下降法，没有本质的差别，在算法实现上也基本⼀样。
因此，本⽂只研究其中⼀种，即FPC算法。
FPC算法虽然实现简单，但其迭代速度慢，效率不⾼。
在此基础上，⽂献[4]做出了改进，提出⼀种Accelerated Proximal Gradient Singular Value Thresholding（APG）算法（该算法是在SVT算法上改进的，本⽂认为FPC和SVT实质上是⼀种算法，故不做区别），能够⼤幅度地提⾼收敛速度。
前⾯提到的⼏种算法，都是从(1-1)中的最⼩化秩的问题出发，经过⼀步步凸松弛得到的。
与上述基本思路不同，⽂献[5]提出了OptSpace算法，它实质上是通过解另⼀种优化问题来实现矩阵完备：min F(W)=∑(i;j)∈Ω∥M ij−W ij∥2s.t.rank(W)=r(1-5)该优化问题旨在找到⼀个秩为r的矩阵，使得该矩阵在对应采样元素的位置上和原始矩阵尽量接近。
在(1-5)中，秩为r的矩阵W可以表⽰成如下分解形式:W=XSY T(1-6)其中，X∈R m×r，Y∈R n×r，S∈R r×r，并且X T X=mI，Y T Y=nI。
注意到(1-6)中的分解与奇异值分解具有类似的形式，X、Y的各列互相正交，但S不是对⾓阵。
如果给定⼀组(X,Y)，则可以找到最佳的S，使得XSY T在对应采样元素的位置上和原始矩阵尽量接近。
因此可以定义如下函数：F(X,Y)=minS∈R r×r∑(i;j)∈Ω∥M ij−(XSY T)ij∥2(1-7)该函数的定义域为D(F)={(X,Y)|X T X=mI,Y T Y=nI,X∈R m×r,Y∈R n×r}。
优化问题(1-5)可以等价为：min F(X,Y)(X,Y)∈D(F)(1-8)函数F的定义域D(F)实质上是两个Grassmann流形的笛卡尔积。
对于在Grass-mann流形上的优化问题，⽂献[6]做了相关研究，并提出了有效的求解算法，⽐如⽜顿法和共轭梯度算法。
⽂献[5]提出的OptSpace算法本质上是求解(1-8)中的优化问题。
仿真结果证明，OptSpace算法相对于FPC、APG算法，能够以更⾼的精度恢复出原始矩阵。
从(1-5)-(1-8)式中可以发现，OptSpace算法似乎要求已知对待恢复矩阵的秩。
但⽂献[7]提出了⼀种从不完全采样矩阵中估计矩阵的秩的算法。
因此，在本⽂中可以假设秩r已知。
1.2本文主要内容本⽂主要对FPC、APG和OptSpace三种算法进⾏了研究，总结了三种算法的基本流程，并通过仿真对三种算法的性能进⾏⽐较，主要⽐较指标为迭代次数、运算时间和恢复精度。
本⽂原创性的⼯作在于，⽂献[5]提出的OptSpace算法，只给出了其⼀般流程，其中包含了两个⼦优化问题，原⽂并没有给出具体算法。
本⽂给出了求解这两个⼦问题的具体算法，完善了OptSpace算法流程。
2FPC算法FPC算法实质上是求解(1-4)中优化问题的⼀种梯度下降求法。
FPC算法的基本思路是：固定某⼀点，在该点对⽬标函数进⾏⼀阶近似，求得最⼩值点，作为下⼀步迭代的固定点，通过这样反复迭代得到全局最⼩值。
记⽬标函数Q(W)=f(W)+µ∥W∥∗f(W)=12∥A(W−b)∥22(2-9)对⽬标函数Q(W)在Z点进⾏⼀阶展开得Qτ(W,Z)=f(Z)+<∇f(Z),W−Z>+τ2∥W−Z∥2F+µ∥W∥∗=τ2∥W−G∥2F+µ∥W∥∗−12τ∥∇f(Z)∥2F(2-10)其中，τ为确定的常数，G的表达式如下：G=Z−τ−1∇f(Z)=Z−τ−1A∗(A(Z)−b)(2-11)在式(2-11)中，A∗表⽰映射A的伴随算⼦。
注意到，映射A实际上是将矩阵中采样位置上的元素提出来排成⼀个向量，⽽A∗是将向量还原到矩阵对应的采样位置上，其余位置上的元素置为0。
由于(2-10)是⼀个关于W凸函数，对于这样⼀个凸优化问题，它的最优解可以表⽰成：Sµτ(G)=U Diag((σ−µ/τ)+)V T(2-12)其中，G的奇异值分解为G=U Diag(σ)V T，(·)+=max{·,0}。
根据上述FPC算法的基本思路以及式(2-12)，算法1给出了FPC算法的⼀般步骤。
算法1FPC算法给定参数µ>0，选取初值W0∈R m×n，对于k=0,1,2,···，按照下述步骤进⾏迭代1:计算G k=X k−(τk)−1A∗(A(X k)−b)2:对G k进⾏奇异值分解，并计算Sµτk(G k)=U Diag((σ−µ/τk)+)V T3:令W k+1=Sµτk(G k)，返回第1步在上述算法中，τk 可以取函数f 的Lipschitz 常数L f 。
对于(2-9)中的函数f (W )=12∥A (W −b )∥22，L f =1。
因此，在上述算法中，可取τk =1。
迭代停⽌的条件为∥W k +1−W k ∥F <ϵ，ϵ为给定精度。
3APG 算法FPC 算法思路简单，易于实现，但是收敛速度⽐较慢。
⽂献[4]对FPC 算法进⾏了加速改进，通过改进下降步长，能够提⾼算法速度。
APG 算法的基本思路与FPC 算法基本⼀致，唯⼀的不同是改进了下降步长。
算法2给出了APG 算法的⼀般步骤。
算法2APG 算法给定参数µ>0，选取初值W 0=W −1∈R m ×n ，初始化步长参数t 0=t −1=0，对于k =0,1,2,···，按照下述步骤进⾏迭代1:令Z k =W k +t k −1−1t k(W k −W K −1)2:计算G k =Z k −(τk )−1A ∗(A (Z k )−b )3:对G k 进⾏奇异值分解，并计算S µτk (G k )=U Diag ((σ−µ/τk )+)V T 4:令W k +1=S µτk (G k )5:计算tk +1=1+√1+4(t k )22，返回第1步同样，在算法2中，取τk =1。
迭代停⽌的条件为∥W k +1−W k ∥F <ϵ，ϵ为给定精度。
4OptSpace 算法⽂献[5]中提出的OptSpace 算法，通过求解下⾯的优化问题来实现矩阵完备：min F (X,Y )=min S ∈Rr ×r∑(i ;j )∈Ω∥M ij −(XSY T )ij ∥2(4-13)由于优化的⽬标函数F (X,Y )的定义域，D (F )={(X ,Y )|X T X =mI,Y T Y =nI,X ∈R m ×r ,Y ∈R n ×r }，为Grassmann 流形，采⽤流形上的梯度下降算法，能够得到该优化问题的解。
注意到，函数F (X,Y )本⾝就是⼀个优化后的最⼩值，没有明确的解析解，但只要给定⼀组(X,Y )，就能求得最⼩化∑(i ;j )∈Ω∥M ij −(XSY T )ij ∥2的S ，进⽽得到F (X,Y )的值。
为后⾯表⽰⽅便，记x =(X,Y )。
4OPTSPACE算法6在给出算法之前，先给出函数F(x)在流形上梯度的表达式。
⾸先，定义⼀个投影算⼦PΩ，∀W∈R m×n，PΩ(W)ij={M ij if(i,j)∈Ω0otherwise(4-14)函数F(x)对于X、Y两部分的梯度分别是grad F(x)X=PΩ(W)(XSY T−M)Y S T−XQ X(4-15)grad F(x)Y=PΩ(W)(XSY T−M)T XS−Y Q Y(4-16)其中，S为给定(X,Y)后最⼩化∑(i;j)∈Ω∥M ij−(XSY T)ij∥2的解，Q X,Q Y∈R r×r分别为：Q X=1mX T PΩ(W)(XSY T−M)Y S T(4-17)Q X=1nY T PΩ(W)(XSY T−M)T XS(4-18)算法3给出OptSpace算法的⼀般流程：算法3OptSpace算法选取初值x0=(X0,Y0)∈D(F)，对于k=0,1,2,···，按照下述步骤进⾏迭代1:针对给定的x k，计算相应的S k2:计算负梯度w k=−grad F(x)3:令x k(t)=x k+t w k(t>0)，求得使F(x k(t))最⼩的步长t4:令x k+1=x k+t w k，返回第⼀步迭代停⽌的条件为∥X k+1S k+1Y Tk+1−X k S k Y T k∥F<ϵ，ϵ为给定精度。