# A New Algorithm for Positive Semidefinite Matrix Completion

Authors: Fangfang Xu; Peng Pan
Journal: Journal of Applied Mathematics
 (2016)

---

## Abstract

Positive semidefinite matrix completion (PSDMC) aims to recover positive semidefinite and low-rank matrices from a subset of entries of a matrix. It is widely applicable in many fields, such as statistic analysis and system control. This task can be conducted by solving the nuclear norm regularized linear least squares model with positive semidefinite constraints. We apply the widely used alternating direction method of multipliers to solve the model and get a novel algorithm. The applicability and efficiency of the new algorithm are demonstrated in numerical experiments. Recovery results show that our algorithm is helpful.


---

## Body

## 1. Introduction

Matrix completion (MC) is the process of recovering the unknown or missing elements of a matrix. Under certain assumptions on the matrix, for example, low-rank or approximately low-rank, the incomplete matrix can be reconstructed very well [1, 2]. Matrix completion is widely applicable in many fields, such as machine learning, statistic analysis, system control, and image and video processing [3], where matrices with low-rank or approximately low-rank are widely used in the model construction.Recently, there have been extensive research on the problems of low-rank matrix completion (LRMC). The affine rank minimization problem consists of finding a matrix of minimum rank that satisfies a given system of linear constraints. However, it is NP-hard (nondeterministic polynomial-time hard) due to the combinatorial nature of the rank function. References [1, 2] showed that the solution of LRMC could be found by solving a nuclear norm minimization problem under some reasonable conditions. The singular value thresholding (SVT) method [4] and fixed point continuation method using approximate singular value decomposition (FPCA) [5] are two well-known algorithms because of their good recoverability, fast speed, and robustness. SVT applied the linearized Bregman iterations to solve the unconstrained nuclear norm regularized linear least squares problem. FPCA used iterations based on an iterative shrinkage-thresholding algorithm and used the continuation technique together with an approximate singular value decomposition procedure to accelerate the algorithm. Reference [6] proposed an accelerated proximal gradient singular value thresholding algorithm. A completely different model was developed in LMaFit [7], which was a nonlinear successive overrelaxation algorithm that only requires solving a linear least squares problem per iteration. More details on LRMC can be found in [1, 8–10] and references therein.In practice, the completed matrix is often required to be positive semidefinite. For example, covariance matrix and its inverse: precision matrix of statistic analysis, are both positive semidefinite. Recently, there have been extensive research on high-dimensional covariance matrix estimation. They all motivate the development of positive semidefinite matrix completion (PSDMC). Reference [11] accomplished the matrix completion task in some special conditions and used the alternating direction method of multipliers (ADMM) [12–15] to solve the model. References [16, 17] proposed new models for nonnegative matrix completion and also used ADMM to solve them.Our main contribution in this work is the development of an efficient algorithm for PSDMC. First of all, we present the nuclear norm regularized linear least squares model with positive semidefinite constraints. Because of its robustness, we choose it as the model of PSDMC in this paper. The structure of the model suggests an alternating minimization scheme, which is very suitable for solving large-scale problems. We give an exact ADMM-based algorithm, whose subproblems are solved exactly. We test the new ADMM-based algorithm on two kinds of problems: random matrix completion problems and random low-rank approximation problems. Numerical experiments show that all our proposed algorithm outputs have satisfactory results. The paper is organized as follows. Section2 presents models and algorithms of PSDMC. Some numerical results are given in Section 3.The following notations will be used throughout this paper. Uppercase (lowercase) letters are used for matrices (column vectors). All vectors are column vectors; the subscript(·)T denotes matrix and vector transposition. Diag(x) denotes a diagonal matrix with x on its main diagonal. 0 is a matrix of all zeros of proper dimension; In stands for the n×n identity matrix. The trace of X∈Rm×n, that is, the sum of the diagonal elements of X, is denoted by tr(X). The Frobenius norm of X∈Rm×n is defined as XF=∑i,j‍|Xi,j|2. The Euclidean inner product between two matrices X∈Rm×n and Z∈Rm×n is defined as X,Z=∑i,j‍(Xi,jZi,j)=tr(X⊤Z). The inequality X⪰0 means that X is semidefinite positive. The equality X=Z means that Xij=Zij for all entries (i,j).


## 2. ADMM-Based Methods for PSDMC

### 2.1. The Model of PSDMC

The matrix completion problem of recovering a positive semidefinite low-rank matrix from a subset of its entries is(1)minrank⁡Xs.t.Xij=Mij,∀i,j∈ΩX⪰0, where X∈Rm×n is the decision variable and Ω is the index set of p known elements of X.LetP be the projection onto the subspace of sparse matrices with nonzeros restricted to the index set Ω; that is,(2)PΩXij=Xij,ifi,j∈Ω,0,otherwise.From the definition ofPΩ, we can reformulate the equality constraint in model (1) in terms of PΩ(X)=PΩ(M). Due to the combinational property of the objective function rank(·), model (1) is NP-hard in general. Inspired by the success of matrix completion under nuclear norm in [1, 2, 10], we use the nuclear norm as an approximation to rank(X) to estimate the optimal solution X∗ of model (1) from the following model: (3)minX∗s.t.PΩX=PΩMX⪰0, where the nuclear norm X∗ of X is defined as the summation of the singular values of X; that is,(4)X∗=∑i=1minm,nσiX, where σi(X) is the ith largest singular value. Moreover, for a positive semidefinite matrix, X∗=I,X=tr(X).If the known elements of the matrixX are noise-free, that is to say, PΩ(M) is reliable, we will directly solve model (3) to conduct SDPMC. On the contrary, if the vector of known elements PΩ(M) is contaminated by noise, the constraints PΩ(X)=PΩ(M) must be relaxed, resulting in the following problem: (5)minXtrX,s.t.PΩX-PΩMF⩽δ,X⪰0, or the nuclear norm regularized linear least squares model with positive semidefinite constraints: (6)minXμtrX+12PΩX-PΩMF2,s.t.X⪰0. Here, δ and μ are given parameters, whose values should be set according to the noise level. When the values of μ and δ are set properly, (5) and (6) are equivalent. Model (6) is usually preferred over (5) for the case of noisy observations. Our algorithms can be extended to treat (5) with minor modifications.Actually, model (6) is especially useful in practice. The reason is that the known information is usually gotten from large surveys and contaminated by sampling error inevitably. In this paper, we choose model (6) as the model to conduct NMC.


### 2.2. An ADMM-Based Method for Model (6)

In this subsection, we present an algorithm developed for model (6). To facilitate an efficient use of ADMM, we introduce one new matrix (splitting) variable Y and consider an equivalent form of model (6): (7)minX,YμtrX+12PΩY-PΩMF2,s.t.X=Y,X⪰0, where X, Y∈Rn×n. The augmented Lagrangian function of model (7) is (8)LX,Y,Π=μtrX+12PΩY-PΩMF2+Π,X-Y+ρ2X-YF2, where Π∈Rn×n is a Lagrangian multiplier. ρ>0 is a penalty parameter.The alternating direction method of multipliers for model (7) is derived by successively minimizing L with respect to X and Y in an alternating fashion; namely,(9a)Xk+1≔argminX⪰0⁡LX,Yk,Πk,(9b)Yk+1≔argminY⁡LXk+1,Y,Πk,(9c)Πk+1≔Πk+γρXk+1-Yk+1,where γ∈(0,1.618).By rearranging the terms of (9a), it is equivalent to(10)minX⪰0⁡μtrX+ρ2X-GkF2,where Gk=Yk-1/ρΠk. Let its Eigenvalue Value Decomposition (EVD) be Gk=UkDiag(wk)UkT, where Uk∈Rn×n and wk∈Rn. Define the shrinkage operator Sν(·) as(11)Sνw=w¯,withw¯i=wi-ν,ifwi-ν>0,0,otherwise.Then(12)Xk+1≔UkDiagSμ/ρwkUkTis an optimal solution of problem (9a).By rearranging the terms of (9b), it is equivalent to (13)minY∈Rn×n12PΩY-PΩMF2+ρ2Y-Xk+1ρΠkF2. Model (13) can be split into two subproblems: (14)min12PΩY-PΩMF2+ρ2PΩY-PΩXk+1ρΠkF2,minPΩ^Y-PΩ^Xk+1ρΠkF2,where Ω^ is the complement of Ω. Finally, we can get the solution of (13) by solving the above two subproblems. (15)Yk+1Ω=1ρ+1PΩM+ρEk,Yk+1Ω^=PΩ^Ek,where Ek=Xk+1/ρΠk.In short, ADMM applied to model (7) yields the iteration:(16a)Xk+1≔UkDiagSμ/ρwkUkT,(16b)Yk+1Ω≔1ρ+1PΩM+ρEk,Yk+1Ω^≔PΩ^Ek,(16c)Πk+1≔Πk+γρXk+1-Yk+1.From the above considerations, we arrive at Algorithm1.Algorithm 1:An exact ADMM-based algorithm of PSDMC.
(1) Input PΩ(M), Im≥0, and tol≥0. (2) Set μ, γ, and ρ≥0. Set X0 and Y0 as random matrices, and Π0 as zero matrices of appropriate sizes. (3) while  notconverge  do (4)  update Xk,Yk,Πk by the formulas (16a), (16b) and (16c).The convergence of ADMM and its variants for convex problems has been studied extensively. The interested reader is referred to [12–15, 18–20], and the references therein. Since model (6) is convex, the convergence of Algorithm 1 to a global optimal solution is guaranteed.


## 2.1. The Model of PSDMC

The matrix completion problem of recovering a positive semidefinite low-rank matrix from a subset of its entries is(1)minrank⁡Xs.t.Xij=Mij,∀i,j∈ΩX⪰0, where X∈Rm×n is the decision variable and Ω is the index set of p known elements of X.LetP be the projection onto the subspace of sparse matrices with nonzeros restricted to the index set Ω; that is,(2)PΩXij=Xij,ifi,j∈Ω,0,otherwise.From the definition ofPΩ, we can reformulate the equality constraint in model (1) in terms of PΩ(X)=PΩ(M). Due to the combinational property of the objective function rank(·), model (1) is NP-hard in general. Inspired by the success of matrix completion under nuclear norm in [1, 2, 10], we use the nuclear norm as an approximation to rank(X) to estimate the optimal solution X∗ of model (1) from the following model: (3)minX∗s.t.PΩX=PΩMX⪰0, where the nuclear norm X∗ of X is defined as the summation of the singular values of X; that is,(4)X∗=∑i=1minm,nσiX, where σi(X) is the ith largest singular value. Moreover, for a positive semidefinite matrix, X∗=I,X=tr(X).If the known elements of the matrixX are noise-free, that is to say, PΩ(M) is reliable, we will directly solve model (3) to conduct SDPMC. On the contrary, if the vector of known elements PΩ(M) is contaminated by noise, the constraints PΩ(X)=PΩ(M) must be relaxed, resulting in the following problem: (5)minXtrX,s.t.PΩX-PΩMF⩽δ,X⪰0, or the nuclear norm regularized linear least squares model with positive semidefinite constraints: (6)minXμtrX+12PΩX-PΩMF2,s.t.X⪰0. Here, δ and μ are given parameters, whose values should be set according to the noise level. When the values of μ and δ are set properly, (5) and (6) are equivalent. Model (6) is usually preferred over (5) for the case of noisy observations. Our algorithms can be extended to treat (5) with minor modifications.Actually, model (6) is especially useful in practice. The reason is that the known information is usually gotten from large surveys and contaminated by sampling error inevitably. In this paper, we choose model (6) as the model to conduct NMC.


## 2.2. An ADMM-Based Method for Model (6)

In this subsection, we present an algorithm developed for model (6). To facilitate an efficient use of ADMM, we introduce one new matrix (splitting) variable Y and consider an equivalent form of model (6): (7)minX,YμtrX+12PΩY-PΩMF2,s.t.X=Y,X⪰0, where X, Y∈Rn×n. The augmented Lagrangian function of model (7) is (8)LX,Y,Π=μtrX+12PΩY-PΩMF2+Π,X-Y+ρ2X-YF2, where Π∈Rn×n is a Lagrangian multiplier. ρ>0 is a penalty parameter.The alternating direction method of multipliers for model (7) is derived by successively minimizing L with respect to X and Y in an alternating fashion; namely,(9a)Xk+1≔argminX⪰0⁡LX,Yk,Πk,(9b)Yk+1≔argminY⁡LXk+1,Y,Πk,(9c)Πk+1≔Πk+γρXk+1-Yk+1,where γ∈(0,1.618).By rearranging the terms of (9a), it is equivalent to(10)minX⪰0⁡μtrX+ρ2X-GkF2,where Gk=Yk-1/ρΠk. Let its Eigenvalue Value Decomposition (EVD) be Gk=UkDiag(wk)UkT, where Uk∈Rn×n and wk∈Rn. Define the shrinkage operator Sν(·) as(11)Sνw=w¯,withw¯i=wi-ν,ifwi-ν>0,0,otherwise.Then(12)Xk+1≔UkDiagSμ/ρwkUkTis an optimal solution of problem (9a).By rearranging the terms of (9b), it is equivalent to (13)minY∈Rn×n12PΩY-PΩMF2+ρ2Y-Xk+1ρΠkF2. Model (13) can be split into two subproblems: (14)min12PΩY-PΩMF2+ρ2PΩY-PΩXk+1ρΠkF2,minPΩ^Y-PΩ^Xk+1ρΠkF2,where Ω^ is the complement of Ω. Finally, we can get the solution of (13) by solving the above two subproblems. (15)Yk+1Ω=1ρ+1PΩM+ρEk,Yk+1Ω^=PΩ^Ek,where Ek=Xk+1/ρΠk.In short, ADMM applied to model (7) yields the iteration:(16a)Xk+1≔UkDiagSμ/ρwkUkT,(16b)Yk+1Ω≔1ρ+1PΩM+ρEk,Yk+1Ω^≔PΩ^Ek,(16c)Πk+1≔Πk+γρXk+1-Yk+1.From the above considerations, we arrive at Algorithm1.Algorithm 1:An exact ADMM-based algorithm of PSDMC.
(1) Input PΩ(M), Im≥0, and tol≥0. (2) Set μ, γ, and ρ≥0. Set X0 and Y0 as random matrices, and Π0 as zero matrices of appropriate sizes. (3) while  notconverge  do (4)  update Xk,Yk,Πk by the formulas (16a), (16b) and (16c).The convergence of ADMM and its variants for convex problems has been studied extensively. The interested reader is referred to [12–15, 18–20], and the references therein. Since model (6) is convex, the convergence of Algorithm 1 to a global optimal solution is guaranteed.


## 3. Numerical Results

In this section, we report on the application of our proposed ADMM-based algorithm to a series of matrix problems to demonstrate its ability. To illustrate the performance of our algorithmic approaches combined with different procedures, we test the following two solvers.(1)
ADMM-VI; Algorithm2 in [11].(2)
ADMM-SDP; Algorithm1, the exact ADMM-based method for model (6).We implement our algorithms in MATLAB. All the experiments are performed on a 2.20 GHz Intel Pentium PC with 6.0 GHz of memory and MATLAB 2012b.


### 3.1. Implementation and Parameters of the Two Solvers

We test the above two solvers on random positive semidefinite matrix problems. We do numerical experiments by the following procedure. Firstly, we create a low-rank positive semidefinite matrixM∈S+n. Secondly, we select a subset of p elements uniformly at random from the n×n elements of M and denote their index set as Ω. The index set of unknown elements is denoted as Ω^. The ratio p/n2 between the number of measurements and the number of entries in the matrix is denoted as “SR” (sampling ratio). We use PΩ(M) to complete X; the result is denoted as X^. Finally, we will check the differences between X^ and its actual value M. (17)err=X^-MF.The most important algorithmic parameters in Algorithm1 are μ, γ, tol, ρ, and the maximal number of iterations Im. In our implementation, we set γ=1.618, tol=2×10-4, ρ=2.5/n, and Im=300. However, the value of μ is very difficult to set since it can be neither too large nor too small. μ is usually chosen to be a moderate value. Three sets of results are computed when μ is equal to three different values. They are shown in Table 1.Table 1
Numerical results of different values ofμ (n=2000, r=20, σ=0.01).
SR μ = 1 0 - 2 μ = 1 0 - 1 μ = 1 μ = 10 Iter Rank Err Iter Rank Err Iter Rank Err Iter Rank Err 0.25 67 297 6.54 e - 03 57 20 2.93 e - 03 44 20 7.79 e - 03 72 20 8.80 e - 02 0.35 50 274 5.61 e - 03 38 20 2.37 e - 03 32 20 7.17 e - 03 51 20 5.91 e - 02 0.45 39 194 3.90 e - 03 28 20 1.97 e - 03 26 20 5.05 e - 03 39 20 4.62 e - 02 0.55 33 137 2.91 e - 03 22 20 1.77 e - 03 21 20 3.53 e - 03 31 20 3.69 e - 02 0.65 29 97 2.26 e - 03 17 20 1.58 e - 03 17 20 3.15 e - 03 26 20 2.97 e - 02 0.75 26 55 1.71 e - 03 15 20 1.43 e - 03 14 20 2.87 e - 03 22 20 2.62 e - 02 0.85 23 20 1.46 e - 03 12 20 1.37 e - 03 11 20 2.75 e - 03 19 20 2.25 e - 02We can use continuation technique employed in [5, 6] to accelerate the convergence of ADMM-SDP. If model (6) is to be solved with the target parameter value μ=μ¯, we propose solving a sequence of model (6) by a decreasing sequence {μ0,μ1,…}. When a new problem, associated with μj+1, is to be solved, the approximate solution for the current problem with μ=μj is used as the starting point. In our numerical experiments in Section 3.2, we set the initial μ0=1 and update μk=μk-1/1.01 at iteration k. A stopping criterion of Algorithm 1 is met as long as all the following three conditions are satisfied:(18a)Xk-Xk-1F≤tol,(18b)Xk-YkF≤tol,(18c)PΩXk-PΩMF≤tol.


### 3.2. Experiments on Random Matrix Completion Problems

The matrixM∈S+m×n with rank r in this subsection is created randomly by the following procedure (see also [5, 7]): two random matrices ML∈Rn×r and MR∈Rn×r with i.i.d. standard Gaussian entries are first generated, Σ∈Rn×n is the noise matrix, and σ>0. Then M=MLMRT+σΣ is assembled. From the creation process of M, its rank is no larger than r.The computational results of positive semidefinite matrix completion are presented in Table2. We can observe that our proposed solver performs as well as ADMM-VI. Both of them can output satisfactory results.Table 2
Numerical results on medium randomly created matrix completion problems.
(n, r, SR) ADMM-VI ADMM-SDP Iter Rank Err Iter Rank Err (500, 10, 0.45) 24 10 2.43 e - 03 29 10 2.62 e - 03 (500, 10, 0.65) 16 10 2.12 e - 03 19 10 2.09 e - 03 (500, 10, 0.85) 10 10 1.90 e - 03 13 10 1.88 e - 03 (2000, 10, 0.45) 27 10 1.86 e - 03 26 10 1.53 e - 03 (2000, 10, 0.65) 18 10 1.19 e - 03 17 10 1.17 e - 03 (2000, 10, 0.85) 11 10 1.24 e - 03 13 10 1.04 e - 03 (2000, 20, 0.45) 27 20 2.12 e - 03 28 20 2.20 e - 03 (2000, 20, 0.65) 16 20 1.84 e - 03 18 20 1.69 e - 03 (2000, 20, 0.85) 11 20 1.46 e - 03 13 20 1.38 e - 03


## 3.1. Implementation and Parameters of the Two Solvers

We test the above two solvers on random positive semidefinite matrix problems. We do numerical experiments by the following procedure. Firstly, we create a low-rank positive semidefinite matrixM∈S+n. Secondly, we select a subset of p elements uniformly at random from the n×n elements of M and denote their index set as Ω. The index set of unknown elements is denoted as Ω^. The ratio p/n2 between the number of measurements and the number of entries in the matrix is denoted as “SR” (sampling ratio). We use PΩ(M) to complete X; the result is denoted as X^. Finally, we will check the differences between X^ and its actual value M. (17)err=X^-MF.The most important algorithmic parameters in Algorithm1 are μ, γ, tol, ρ, and the maximal number of iterations Im. In our implementation, we set γ=1.618, tol=2×10-4, ρ=2.5/n, and Im=300. However, the value of μ is very difficult to set since it can be neither too large nor too small. μ is usually chosen to be a moderate value. Three sets of results are computed when μ is equal to three different values. They are shown in Table 1.Table 1
Numerical results of different values ofμ (n=2000, r=20, σ=0.01).
SR μ = 1 0 - 2 μ = 1 0 - 1 μ = 1 μ = 10 Iter Rank Err Iter Rank Err Iter Rank Err Iter Rank Err 0.25 67 297 6.54 e - 03 57 20 2.93 e - 03 44 20 7.79 e - 03 72 20 8.80 e - 02 0.35 50 274 5.61 e - 03 38 20 2.37 e - 03 32 20 7.17 e - 03 51 20 5.91 e - 02 0.45 39 194 3.90 e - 03 28 20 1.97 e - 03 26 20 5.05 e - 03 39 20 4.62 e - 02 0.55 33 137 2.91 e - 03 22 20 1.77 e - 03 21 20 3.53 e - 03 31 20 3.69 e - 02 0.65 29 97 2.26 e - 03 17 20 1.58 e - 03 17 20 3.15 e - 03 26 20 2.97 e - 02 0.75 26 55 1.71 e - 03 15 20 1.43 e - 03 14 20 2.87 e - 03 22 20 2.62 e - 02 0.85 23 20 1.46 e - 03 12 20 1.37 e - 03 11 20 2.75 e - 03 19 20 2.25 e - 02We can use continuation technique employed in [5, 6] to accelerate the convergence of ADMM-SDP. If model (6) is to be solved with the target parameter value μ=μ¯, we propose solving a sequence of model (6) by a decreasing sequence {μ0,μ1,…}. When a new problem, associated with μj+1, is to be solved, the approximate solution for the current problem with μ=μj is used as the starting point. In our numerical experiments in Section 3.2, we set the initial μ0=1 and update μk=μk-1/1.01 at iteration k. A stopping criterion of Algorithm 1 is met as long as all the following three conditions are satisfied:(18a)Xk-Xk-1F≤tol,(18b)Xk-YkF≤tol,(18c)PΩXk-PΩMF≤tol.


## 3.2. Experiments on Random Matrix Completion Problems

The matrixM∈S+m×n with rank r in this subsection is created randomly by the following procedure (see also [5, 7]): two random matrices ML∈Rn×r and MR∈Rn×r with i.i.d. standard Gaussian entries are first generated, Σ∈Rn×n is the noise matrix, and σ>0. Then M=MLMRT+σΣ is assembled. From the creation process of M, its rank is no larger than r.The computational results of positive semidefinite matrix completion are presented in Table2. We can observe that our proposed solver performs as well as ADMM-VI. Both of them can output satisfactory results.Table 2
Numerical results on medium randomly created matrix completion problems.
(n, r, SR) ADMM-VI ADMM-SDP Iter Rank Err Iter Rank Err (500, 10, 0.45) 24 10 2.43 e - 03 29 10 2.62 e - 03 (500, 10, 0.65) 16 10 2.12 e - 03 19 10 2.09 e - 03 (500, 10, 0.85) 10 10 1.90 e - 03 13 10 1.88 e - 03 (2000, 10, 0.45) 27 10 1.86 e - 03 26 10 1.53 e - 03 (2000, 10, 0.65) 18 10 1.19 e - 03 17 10 1.17 e - 03 (2000, 10, 0.85) 11 10 1.24 e - 03 13 10 1.04 e - 03 (2000, 20, 0.45) 27 20 2.12 e - 03 28 20 2.20 e - 03 (2000, 20, 0.65) 16 20 1.84 e - 03 18 20 1.69 e - 03 (2000, 20, 0.85) 11 20 1.46 e - 03 13 20 1.38 e - 03



---

*Source: 1659019-2016-10-25.xml*