低秩半正定矩阵最小二乘恢复算法

余宏伟 蒋轶

摘 要： 多年来矩阵恢复一直是学术界的一个热门研究课题，它被广泛应用于多个技术领域，如计算机视觉、图像恢复以及推荐系统等。考虑其中一种特殊且十分重要的矩阵恢复，即半正定矩阵恢复。通过将此类矩阵恢复问题与基于测距的网络定位问题类比，构造了基于最小二乘的优化模型，运用顺序凸规划（sequential convex programming/SCP）算法，可以高效并精确地求解此问题，从而将缺失矩阵较为精准地还原为全矩阵。仿真结果证明相比于目前文献中已有矩阵恢复算法，提出的算法具有更好的恢复性能。

关键词： 低秩矩阵;半正定矩阵;核范数最小化;最小二乘法;顺序凸规划

中图分类号： TN92

文献标志码： A

文章编号：1007-757X（2019）06-0001-03

Abstract： Matrix completion has been an active topic for more than a decade in academic research due to its widely applications in many fields， such as computer vision， image recovery and recommendation system. In this paper， we study a special and important case， low-rank positive semi-definite （PSD） matrix completion. By linking the matrix completion problem to a network localization problem， we propose an optimization model based on the least square criterion. Using the sequential solution programming （SCP） algorithm， we can efficiently reconstruct the PSD matrix with high precision. Simulation results show that the proposed method has more superior performance than the state-of-art algorithms.

Key words： Low-rank matrix; PSD matrix; Nuclear norm minimization; Least square; Sequential solution programming

0 引言

矩阵恢复技术要求能够从部分甚至是受到噪声污染的矩阵元素观测中精确地还原整个矩阵。在具有里程碑意义的论文[1]中，Candes等人证明，如果某矩阵是一个低秩矩阵，且其中被观测到的元素集合满足一定的条件，那么通过一个简单的凸规划，就可以以较高的概率将其恢复。矩阵恢复的应用范围非常广泛，如计算机视觉、图像恢复、定位[2]、推荐系统（Netflix）等等。

Candes等人提出通过最小化核范数的方法来恢复整个矩阵[1、3]。特别地，在满足一定非相干条件的情况下，最小化核范数法可以逼近最优解[4]。[1]中应用半正定规划解算器SDPT3求解核范数最小化问题，该计算正如SeDuMi一样都是基于内点法，当矩阵规模很大时需要求解一个巨大的线性方程组计算牛顿方向，一般来说当矩阵维度大于100时，SDPT就无能为力了[5]。后来的很多学者针对此优化问题的求解进行了大量的工作。如[5]中提出了Singular Value Thresholding（SVT）算法，该算法在迭代求解过程中产生一组矩阵序列{Xk，Yk}，每一步迭代只需Yk对进行SVD分解并对奇异值作阈值化操作，最终{Xk}会收敛逼近核范数最小化问题的解。针对存在观测噪声情形下的半正定矩阵恢复，论文[6]提出了一种基于Alternating Direction Method of Multipliers （ADMM）的算法，引入了一个与待恢复矩阵X相等的矩阵变量Y，每次迭代更新X，Y，更新值都存在闭式表达，求解过程简洁高效。类似的工作还有很多[7、8、9]，但研究重点大都是集中在如何更高效的求解Candes所提出的核范数最小化模型，鲜有矩阵恢复性能上的提升。

本文研究一种极为特殊和重要的矩阵，即半正定矩阵的恢复问题。我们将此类矩阵恢复问题与基于测距的网络定位问题[10]类比。从网络定位的视角出发，本文构造了基于最小二乘的半正定矩阵恢复优化模型，凭借我们提出的SCP算法，该优化问题可以被高效精确求解，从而高精度恢复整个半正定矩阵。

本文符号表示说明如下。粗体小写和大写字母分别表示向量和矩阵。A，A分别表示矩阵A的Frobenius范数和核范数。A⊙B表示矩阵A、B的哈达玛积，即对应元素相乘。vec（X）表示对矩阵X按列向量化。

1 问题描述与经典算法

1.1 问题描述

矩阵恢复问题是将整个矩阵从它的部分观测元素集恢复出来。定义A∈Rm×n为待恢复的矩阵。Ω[m]×[n]代表观测到的矩阵A的元素集合，其中[m]表示集合{1，2，…，m}. W∈Rm×n表示掩模矩阵如式（1）。

1.2 经典算法

矩阵恢复的一种经典思路就是从最小化矩阵的秩角度考虑，此时如果并不考虑噪声，那么就形成了以下优化问题，如式（2）。

众所周知，上式优化问题是一个NP-hard问题。为了求解此问题，Candes将其放松成一个最小化矩阵核范數的凸优化问题[1]，如式（3）。

此優化模型较为容易求解，且能够较为精确的从部分无噪声观测集中恢复出整个矩阵A。而且，此模型也可以实现噪声环境下的矩阵恢复，只需稍作改动[1][2]，如式（4）。

如果待恢复矩阵是半正定矩阵，只需要增加一项半正定约束：A≥0，而A*=trace（A）。有很多学者在此基础上进行了大量的工作，但着力点大都只在于如何更高效地求解式（3）中的优化问题，而鲜有恢复性能上的提升。

本文研究一种在通信等应用中极为重要的矩阵，即半正定矩阵的恢复问题。接下来我们将展示我们的方法较经典核范数方法在恢复性能方面有明显的突破。

2 最小二乘法矩阵恢复

如何求解式（5）也是一个具有挑战性的问题。我们在论文[10]中提出了运用顺序凸规划（SCP）算法求解网络定位问题。优化问题[4]实际就是一个类似的定位问题，同样可以用SCP算法求解。该算法需要一组初始值启动牛顿迭代。此处我们直接利用核范数最小化算法的结果构造初始值，在此基础上，可以显著提高矩阵恢复性能。具体构造方法如下。

将式（6）和（7）带入SCP算法，式（5）可以被精确且高效求解。进而可以重构待恢复的半正定矩阵。为了表述清晰和完整性，我们将整个矩阵恢复算法的具体步骤总结在算法1中。值得注意的是，由于（5）不是一个关于X的凸函数，2f可能存在负特征值。因此，我们将海森矩阵的逆（2f）-1用（2f）-代替，即将前者中的负数特征值用0替换而得到的新矩阵。这是SCP算法的核心技巧。

本节将通过仿真验证上述算法的恢复性能。

此处，A∈R20，rank（A）=4。这些虚拟的节点坐落在20×20×20×20的空间中。观测噪声服从高斯分布N（0，1）。本仿真比较在删除不同对数矩阵元素的情形下，矩阵恢复的性能。注意元素删除都是成对删除的，即如果Aij缺失，那么Aji也一定缺失。对于每一种情形，我们都进行了500次蒙特卡罗仿真，每次从N+（N-1）+…+1=210对元素随机选取指定对数的矩阵元素删除。参考论文[2]中的做法，我们也定义矩阵恢复误差为A^-AFAF，并且如果该误差小于1%就认为这是一次成功的矩阵恢复。

比较了无观测噪声和有观测噪声情况下的半正定矩阵恢复成功率高低，如图1和图2所示。

在有观测噪声情况下，直到60对（约29%）元素删除，两种方法都能完全恢复出原矩阵。当80对（约38%）元素删除时，本文提出的最小二乘算法仍能几乎100%恢复，而核范数最小化算法已经出现了明显的失败率，等到删除100对（约48%）时，该方法成功率只剩约60%，已经不太可靠了。而我们的方法依然保持约95%以上成功率，即使是删除120对（约57%）仍然有约60%成功率，此时，核范数最小化算法已经完全失效了。从图1中可以看出，在无噪声情况下，两种方法恢复性能较有噪声情形性能都有所提升，且最小二乘法依然明显优于核范数最小化算法。本文的方法在删除140对元素时也无法实现矩阵恢复，原因如下：半正定矩阵A的自由度是N+（N-1）+…+（N-r+1）=74，要想恢复A，它的210对矩阵元素至多能删除210-74=136对元素。

比较了在无噪声和有噪声情形下半正定矩阵恢复的平均误差，如图3和图4所示。

从图中可以看出，随着删除的矩阵元素对数增加，矩阵恢复性能逐渐下降，但是可以看出无论有无噪声，最小二乘法的平均恢复误差始终显著低于核范数最小化算法的平均误差。

4 总结

本文主要研究低秩半正定矩阵恢复问题，我们从 “测距”网络定位的视角出发，构造出基于最小二乘的优化模型，然后应用顺序凸规划（SCP）算法，可以精确求解并重构待恢复矩阵。仿真试验表明，该方法无论是在无观测噪声还是存在观测噪声情形下都能较为精确的恢复目标矩阵，在核范数最小化算法的基础上获得性能的显著提升。

参考文献

[1] Candes Emmanuel J， Recht Benjamin. Exact matrix completion via convex optimization[J]. Foundations of Computational mathematics， 2009（6）： 717-773.

[2] Dokmanic Ivan， Parhizkar Reza， Ranieri Juri， et al. Euclidean distance matrices： essential theory， algorithms， and applications[J]. IEEE Signal Processing Magazine， 2015，32（6）： 12-30.

[3] Candes Emmanuel J， Plan Yaniv. Matrix completion with noise[J]. Proceedings of the IEEE， 2010，98（6）： 925-936.

[4] Candes Emmanuel J， Tao Terence. The power of convex relaxation： Near-optimal matrix completion[J]. IEEE Transactions on Information Theory， 2010， 56（5）：2053-2080.

[5] Cai Jian-Feng， Candes Emmanuel J， Shen Zuowei. A singular value thresholding algorithm for matrix completion[J]. SIAM Journal on Optimization， 2010， 20（ 4）： 1956-1982.

[6] Xu， Fangfang and Pan， Peng， "A New Algorithm for Positive Semidefinite Matrix Completion，" Journal of Applied Mathematics， vol. 3， pp. 1-5， 2016.

[7] Ma S， Goldfarb D， Chen L. Fixed point and Bregman iterative methods for matrix rank minimization[J]. Mathematical Programming， 2011， 128（1-2）： 321-353.

[8] Fornasier Massimo， Rauhut Holger， Ward Rachel. Low-rank matrix recovery via iteratively reweighted least squares minimization[J]. SIAM Journal on Optimization， 2011， 21（4）： 1614-1640.

[9] Chen Caihua， He Bingsheng， Yuan Xiaoming. Matrix completion via an alternating direction method[J]. IMA Journal of Numerical Analysis， 2012， 3（1）：227-245.

[10] Yu， Hongwei and Jiang， Yi， "Maximum likelihood network localization using range estimation and GPS measurements，" 2017 9th International Conference on Wireless Communications and Signal Processing （WCSP）， Nanjing， 2017. Oct， pp. 1-6.

（收稿日期： 2018.12.27）