A Relaxed Interior Point Method for Low-Rank Semidefinite Programming Problems with Applications to Matrix Completion

A new relaxed variant of interior point method for low-rank semidefinite programming problems is proposed in this paper. The method is a step outside of the usual interior point framework. In anticipation to converging to a low-rank primal solution, a special nearly low-rank form of all primal iterates is imposed. To accommodate such a (restrictive) structure, the first order optimality conditions have to be relaxed and are therefore approximated by solving an auxiliary least-squares problem. The relaxed interior point framework opens numerous possibilities how primal and dual approximated Newton directions can be computed. In particular, it admits the application of both the first- and the second-order methods in this context. The convergence of the method is established. A prototype implementation is discussed and encouraging preliminary computational results are reported for solving the SDP-reformulation of matrix-completion problems.

1. Introduction.We are concerned with an application of an interior point method (IPM) for solving large, sparse and specially structured positive semidefinite programming problems (SDPs).
Let SR n×n denote the set of real symmetric matrices of order n and let U • V denote the inner product between two matrices, defined by trace(U T V ).Consider the standard semidefinite programming (SDP) problem in its primal form where A i , C ∈ SR n×n and b ∈ R m are given and X ∈ SR n×n is unknown and assume that matrices A i , i = 1, 2, . . ., m are linearly independent, that is m i=1 d i A i = 0 implies d i = 0, i = 1, . . ., m.The dual form of the SDP problem associated with (1.1) is: where y ∈ R m and S ∈ SR n×n .
The number of applications which involve semidefinite programming problems as a modelling tool is already impressive [29,32] and is still growing.Applications include problems arising in engineering, finance, optimal control, power flow, various SDP relaxations of combinatorial optimization problems, matrix completion or other applications originating from modern computational statistics and machine learning.Although the progress in the solution algorithms for SDP over the last two decades was certainly impressive (see the books on the subject [2,13]), the efficient solution of general semidefinite programming problems still remains a computational challenge.
Among various algorithms for solving (linear) SDPs, interior point methods stand out as reliable algorithms which enjoy guaranteed polynomial worst-case complexity and usually provide accurate solutions within reasonable time.However, when sizes of SDP instances grow, traditional IPMs which require computing exact Newton search directions hit their limits.Indeed, the effort required by the linear algebra in (standard) IPMs may grow as fast as O(n 6 ).
Although there exists a number of alternative approaches to interior point methods, such as for example [8,9,23], which can solve certain SDPs very efficiently, they usually come without polynomial worst case complexity guarantees, a feature offered only by IPMs.Therefore there is a need to develop faster IPM-based techniques which could preserve some of the excellent theoretical properties of these methods, but compromise on the other features in quest for practical computational efficiency.Customized IPM methods have been proposed for special classes of problems.They take advantage of sparsity and structure of the problems, see e.g.[4,5,18,24,27,30] and the references in [1].
In this paper we focus on problems in which the primal variable X is expected to be low-rank at optimality.Such situations are common in relaxations of combinatorial optimization problems [5], for example in maximum cut problems [19], as well as in matrix completion problems [11], general trust region problems and quadratically constrained quadratic problems in complex variables [26].We exploit the structure of the sought solution and relax the rigid structure of IPMs for SDP.In particular we propose to weaken the usual connection between the primal and dual problem formulation and exploit any special features of the primal variable X.
Rank plays an important role in semidefinite programming.For example, every polynomial optimization problem has a natural SDP relaxation, and this relaxation is exact when it possesses a rank-1 solution [26].On the other hand, for any general problem of the form (1.1), there exists an equivalent formulation where an additional bound r on the rank of X may be imposed as long as r is not too small [9].More specifically, under suitable assumptions, there exists an optimal solution X * of (1.1) with rank r satisfying r(r + 1)/2 ≤ m.There have been successful attempts to identify low rank submatrices in the SDP matrices and eliminate them with the aim to reduce the rank and hence the difficulty of solving an SDP.A technique called facial reduction [22] has been analysed and demonstrated to work well in practice.Interestingly, when positive semidefinite programs are solved using interior-point algorithms, then because of the nature of logarithmic barrier function promoting the presence of nonzero eigenvalues, the primal variable X typically converges to a maximum-rank solution [20,26].However, in this paper we aim at achieving the opposite.We want to design an interior point method which drives the generated sequence of iterates to converge to a low-rank solution.We assume that constraint matrices are sparse and we search for a solution X of rank r of the form X = U U T with U ∈ R n×r .Special low-rank structure of X may be imposed directly in problem (1.1), but this excludes the use of an interior point algorithm (which requires all iterates X to be strictly positive definite).Burer and Monteiro [8,9] and their followers [6,7] have used such an approach with great success.Namely, they have substituted U U T for X in (1.1) and therefore have replaced it with the following nonlinear programming problem with U ∈ R n×r .Although such transformation removes the difficult positive definiteness constraint (it is implicit as X = U U T ), the difficulty is shifted elsewhere as both the objective and constraints in (1.3) are no longer linear, but instead quadratic and in general non-convex.In comparison with a standard IPM the method proposed in [8,9] and applied to solve large-scale problems enjoys substantially reduced memory requirements and very good efficiency and accuracy.However, due to nonconvexity of (1.3), local methods may not always recover the global optimum.In [6,7] authors showed that, despite the non-convexity, first-and second-order necessary optimality conditions are also sufficient, provided that rank r is large enough and constraints satisfy some regularity conditions.That is, when applied to several classes of SDPs, the low-rank Burer-Monteiro formulation is very unlikely to converge to any spurious local optima.
In this paper we propose a different approach.We would like to preserve as many of the advantageous properties of interior point methods as possible and expect to achieve it by (i) working with the original problem (1.1) and (ii) exploiting the lowrank structure of X. Knowing that at optimality X is low-rank we impose a special form of the primal variable throughout the interior point algorithm with U ∈ R n×r , for a given r > 0 and µ denoting the barrier term.Hence X is full rank (as required by IPM), but approaches the low-rank matrix as µ goes to zero.Imposing such special structure of X offers an advantage to an interior point algorithm: it can work with an object of size nr rather than a full rank X of size n 2 .We have additionally considered an adaptive choice of r assuming that this rank may not be known a priori.Indeed, the method can start with r equal to 1 or 2 and gradually increase r to the necessary minimum rank (target rank).Remarkably, the method can also handle problems with nearly-low-rank solution, as the primal variable is not assumed to be low-rank along the iterations, but it is gradually pushed to a low-rank matrix.Finally, the presence of the perturbation term µI allows to deal with possibly noisy right-hand side b as well.We also further relax the rigid IPM structure.Starting from a dual feasible approximation, we dispose of dual slack variable S and avoid computations which would involve large Kronecker product matrices of dimension n 2 × n 2 (and that in the worst case might require up to O(n 6 ) arithmetic operations).
The paper is organised as follows.After a brief summary of notation used in the paper provided at the end of this section, in Section 2 we present the general framework and deliver some theoretical insights into the proposed method.In Section 3 we explain the mechanism which allows to adaptively reveal the rank of the minimum rank solution matrix X.The proposed approach offers significant flexibility in the way how Newton-like search directions are computed.They originate from a solution of a least squares problem.We see it in detail in Section 4. Next, in Section 5 we discuss the properties of low-rank SDPs arising in matrix completion problems and in Section 6 we present preliminary computational results obtained with a prototype Matlab implementation of the new algorithm.Finally, in Section 7 we give our conclusions.Appendix A contains some notes on the Kronecker product of two matrices and on matrix calculus.
Notation.The norm of the matrix associated with the inner product between two matrices U • V = trace(U T V ) is the Frobenius norm, written U F := (U • U ) 1/2 , while U 2 denotes the L 2 -operator norm of a matrix.Norms of vectors will always be Euclidean.The symbol I p denotes the identity matrix of dimension p × p.
Let A be the linear operator A : SR n → R m defined by Moreover, let A T denote the matrix representation of A T with respect to the standard bases of R n , that is and where mat is the "inverse" operator to vec (i.e., mat(vec 2. Relaxed interior point method for low-rank SDP.Interior point methods for semidefinite programming problems work with the perturbed first-order optimality conditions for problems (1.1)-(1.2) given by: A general IPM involves a triple (X, y, S), performs steps in Newton direction for (2.1), and keeps its subsequent iterates in a neighbourhood of the central path [2,13].The convergence is forced by gradually reducing the barrier term µ.However, having in mind the idea of converging to a low-rank solution, we find such a structure rather restrictive and wish to relax it.This is achieved by removing explicit S from the optimality conditions and imposing a special structure of X. Substituting S = C − A T y from the first equation into the third one, we get Next, following the expectation that at optimality X has rank r, we impose on X the following special structure with U ∈ R n×r , for a given r > 0. We do not have any guarantee that there exists a solution of (2.2) with such a structure, but we can consider the least-square problem: min where F r µ (U, y) : R n×r × R m → R n 2 +m is given by , µ > 0. (2.5) The nonlinear function F r µ (U, y) has been obtained substituting X = µI n + U U T in (2.2) after vectorization of the second block.The associated system F r µ (U, y) = 0 is overdetermined with (m + n 2 ) equations and (nr + m) unknowns (U, y).In the following, for the sake of simplicity, we identify R n×r × R m with R nr+m .
It is worth mentioning at this point that the use of least-squares type solutions to an overdetermined systems arising in interior point methods for SDP was considered in [25,14].Its primary objective was to avoid symmetrization when computing search directions and the least-squares approach was applied to a standard, complete set of perturbed optimality conditions (2.1).
We propose to apply to problem (2.4) a similar framework to that of interior point methods, namely: Fix µ, iterate on a tuple (U, y), and make steps towards a solution to (2.4).This opens numerous possibilities.One could for example compute the search directions for both variables at the same time, or alternate between the steps in U and in y.
Bearing in mind that (2.1) are the optimality conditions for (1.1) and assuming that a rank r optimal solution of (1.1) exists, we will derive an upper bound on the optimal residual of the least-squares problem (2.4).Assume that a solution (X * , y * , S * ) of the KKT conditions exists such that X * = U * (U * ) T , U * ∈ R n×r , that is Then evaluating (2.5) at (U * , y * ) and using (2.6) we get Consequently, we obtain the following upper bound for the residual of the least-squares problem (2.4): where Assuming to have an estimate of ω * we are now ready to sketch in Algorithm 1 the general framework of a new relaxed interior point method.
To start the procedure we need an initial guess (U 0 , y 0 ) such that U 0 is full column rank and S 0 = C − A T y 0 is positive definite, and an initial barrier parameter µ 0 > 0. At a generic iteration k, given the current barrier parameter µ k > 0, we compute an approximate solution (U k , ȳk ) of (2.4) such that φ µ k (U k , ȳk ) is below µ 2 k ω * .Then, the dual variable y k and the dual slack variable S k are updated as follows: with α k ∈ (0, 1] such that S k remains positive definite.We draw the reader's attention to the fact that although the dual variable S does not explicitly appear in optimality conditions (2.2) or (2.5), we do maintain it as the algorithm progresses and make sure that (S k , y k ) remains dual feasible.Finally, to complete the major step of the algorithm, the barrier parameter is reduced and a new iteration is performed.
Note that so far we have assumed that there exists a solution to (1.1) of rank r.In case such a solution does not exist the optimal residual of the least-squares problem is not guaranteed to decrease as fast as µ 2 k .This apparently adverse case can be exploited to design an adaptive procedure that increases/decreases r without requiring the knowledge of the solution's rank.This approach will be described in Section 3.
Algorithm 1 General framework of the Relaxed Interior Point algorithm for solving low-rank SDP input: (2.9) 3:

5:
Set µ k = σµ k−1 6: end for In the remaining part of this section we state some of the properties of the Algorithm which are essential to make it work in practice.
First we note that dual constraint is always satisfied by construction and the backtracking process at Line 3 is well-defined.This is proved in Lemma 4 of [4] which is repeated below for sake of reader's convenience.
Lemma 2.1.Let ∆S be computed in Step 4 of Algorithm 1 at iteration k and S k−1 be computed at the previous iteration k − 1.Then, there exists k−1 , the thesis holds with .
Note that if backtracking is needed (i.e.α k < 1) to maintain the positive definiteness of the dual variable, then after updating S k in Step 5 the centrality term may increase and it is not guaranteed to remain below γµ k .Indeed, by setting (2.10) that is the centrality measure may actually increase along the iterations whenever α k does not approach one as µ k goes to zero.In the following we analyse the convergence properties of Algorithm 1 when this adverse situation does not occur, namely under the following assumption: To the best of authors knowledge, it does not seem possible to demonstrate that eventually α k approaches one.This is because we impose a special form of X in (2.3) and make only a weak requirement regarding the proximity of the iterate to the central path: with γ possibly greater than one.However, we will provide computational evidence that this seemingly strong Assumption 1 does hold in the numerical tests.Proposition 2.2.Let Assumption 1 hold.Assume that a solution of rank r of problem (1.1) exists and the sequence {U k , y k } generated by Algorithm 1 is bounded.Then, any limit point Proof.Assume for the sake of simplicity that the whole sequence is converging to which implies X † S † = 0 and by construction ensures that S † is positive semidefinite being a limit point of a sequence of positive definite matrices.
From the previous proposition it follows that (X † , y † , S † ) solves (2.1).Moreover, X † has rank r, unless U † is not full column rank.This situation can happen only in the case (1.1) admits a solution of rank smaller than r.In what follows for sake of simplicity we assume that the limit point U † is full column rank.
Remark 2.3.It is worth observing that due to the imposed structure of matrices (2.3) all iterates X k are full rank, but asymptotically they approach rank r matrix.Moreover, the minimum distance of X k to a rank r matrix is given by µ k , i.e., and the primal infeasibility is bounded by γµ k .This allows us to use the proposed methodology also when the sought solution is close to a rank r matrix ("nearly lowrank") and/or some entries in vector b are corrupted with a small amount of noise.
3. Rank updating/downdating.The analysis carried out in the previous section requires the knowledge of γ ≥ √ ω * and of the rank r of the sought solution.As the scalar γ is generally not known, at a generic iteration k the optimization method used to compute an approximate minimizer of (2.4) is stopped when a chosen first-order criticality measure ψ µ goes below the threshold η 2 µ k where η 2 is a strictly positive constant.This way, the accuracy in the solution of (2.4) increases as µ k decreases, but the magnitude of the optimal residual is not checked.We considered Regarding the choice of the rank r, there are situations where the rank of the sought solution is not known.Below we describe a modification of Algorithm 1 where, starting from a small rank r, the procedure adaptively increases/decreases it.This modification is based on the observation that if a solution of rank r exists the iterative procedure used in Step 2, should provide a sequence {U k } such that the primal infeasibility also decreases with µ k .Then, at each iteration the ratio is checked.If this ratio is larger than η 1 , where η 1 is a given constant in (σ, 1) and σ is the constant used to reduce µ k , then the rank r is increased as the procedure has not been able to provide the expected decrease in the primal infeasibility.After an update of rank, the parameter µ k is not changed and one column is appended to the current U k .As a safeguard, also a downdating strategy can be implemented.In fact, if after an increase of rank, we still have ρ k > η 1 then we come back to the previous rank and inhibit rank updates in all subsequent iterations.This is detailed in Algorithm 2 where we borrowed the Matlab notation.Variable update r is an indicator specifying if at the previous iteration the rank was increased (update r = up), decreased (update r = down) or left unchanged (update r = unch).
The initial rank r should be chosen as the rank of the solution (if known) or as a small value (say 2 or 3) if it is unknown.The dimension of the initial variable U 0 is then defined accordingly.Since given ǫ and σ the number of iterations to satisfy µ k < ǫ at Line 6 is predefined, the number of rank updates is predefined as well.Therefore, if an estimate of the solution rank is known, one should use it in order to define a suitable initial r.
Algorithm 2 Relaxed Interior Point Algorithm for Low Rank SDP (IPLR) input: The initial rank r, initial (y 0 , U 0 ) with U 0 ∈ R n×r and y 0 ∈ R m such that Find an approximate minimizer (U k , ȳk ) of φ µ k (U, y) such that Compute ρ k given in (3.1). 10: if update r = unch then 12: Set r = r + 1 and update r = up 13: 14: else if update r = up then 15: Set r = r − 1 and update r = down 16: end if end if 25: end for 4. Solving the nonlinear least-squares problem.In this section we investigate the numerical solution of the nonlinear least-squares problem (2.4).
Following the derivation rules recalled in Appendix A, we compute the Jacobian matrix J µ k ∈ R (n 2 +m)×(nr+m) of F r µ k which takes the following form: where and Π nr ∈ R nr×nr is the unique permutation matrix such that vec(B T ) = Π nr vec(B) for any B ∈ R n×r , see Appendix A.
In order to apply an iterative method for approximately solving (2.4) we need to perform the action of J T µ k on a vector to compute the gradient of φ µ k .The action of J µ k on a vector is also required in case one wants to apply a Gauss-Newton approach (see Section 4.3).In the next section we will discuss how these computations are carried out.
4.1.Matrix-vector products with blocks of J µ k .First, let us denote the Jacobian matrix blocks as follows: 2) Below we will show that despite J µ k blocks contain matrices of dimension n 2 × n 2 , matrix-vector products can be carried out without involving such matrices and the sparsity of the constraint matrices can be exploited.We will make use of the properties of the Kronecker product (A.1)-(A.3)and assume that if v ∈ R nr and z ∈ R n 2 then mat(v) ∈ R n×r and mat(z) ∈ R n×n .
• Let v ∈ R nr and w ∈ R m and let us consider the action of J 11 and J T 11 on v and w, respectively: and (4.7) • Let v ∈ R nr and z ∈ R n 2 and let us consider the action of J T 21 J 21 and J T 21 on v and z, respectively: and • Let w ∈ R m and z ∈ R n 2 and let us consider the action of J 22 and J T 22 on w and z, respectively: and with Z = (µI n + U U T )mat(z).(4.12) 4.2.Computational effort per iteration.The previous analysis shows that we can perform all products involving Jacobian's blocks handling only n × n matrices.Moreover, if matrices A i are indeed very sparse, as is the case in matrix completion problems, their structure can be exploited in the matrix-products in (4.7) and (4.10).Moreover, only few elements of matrices V in (4.6) and Z in (4.12) need to be involved when products (4.5) and (4.11) are computed, respectively.More precisely, denoting with nnz(A) the number of nonzero entries of A, we need to compute nnz(A) entries of V and Z defined in (4.6) and (4.12), respectively.Noting that mat(v) ∈ R n×r and U T ∈ R r×n , the computation of the needed entries of V amounts to (O(nnz(A)r) flops.Regarding Z, the computation of the intermediate matrix Ŵ = U T mat(z) ∈ R r×n costs O(n 2 r) flops and nnz(A) entries of U Ŵ requires O(nnz(A)r) flops.
In table 4.2 we provide the estimate flop counts for computing various matrixvector products with the blocks of Jacobian matrix.We consider the products that are relevant in the computation of the gradient of φ µ k and in handling the linearalgebra phase of the second order method which we will introduce in the next section.From the table, it is evident that the computation of the gradient of φ µ k requires O(max{nnz(A), n 2 }r + m) flops.Then, under the following assumptions: (i) nnz(A) = O(n 2 ), (ii) a first-order method is applied to problem (2.4), (iii) at most O(n) inner iterations of such method are performed, the complexity of each iteration of Algorithm 1 is O(n 3 r + nm), and it is significantly smaller than O(n 6 ) required by a general purpose interior-point solver [2,13] or O(n 4 ) needed by the specialized approach for nuclear norm minimization [27].Apart from all operation listed above the backtracking along ∆S needs to ensure that S k is positive definite (Algorithm 1, Step 4) and this is verified by computing the Cholesky factorization of the matrix S k−1 + α k ∆S, for each trial steplength α k .If the dual matrix is sparse, i.e. when matrices A i , i = 1, . . ., m and C share the sparsity patterns [31], a sparse Cholesky factor is expected.Note that the structure of dual matrix does not change during the iterations, hence reordering of S 0 can be carried out once at the very start of Algorithm 2 and then may be reused to compute the Cholesky factorization of S k−1 + α k ∆S at each iteration.

Nonlinear
Gauss-Seidel approach.The crucial step of our interior point framework is the computation of an approximate solution of the nonlinear leastsquares problem (2.4).To accomplish the goal, a first order approach as well as a Gauss-Newton method can be used.However, in this latter case the linear algebra phase becomes an issue, due to the large dimension of the Jacobian.Here, we propose a Nonlinear Gauss-Seidel method and show how to handle the linear algebra phase.In this approach, at Step 3 of Algorithm 2 we adopt the process detailed in Algorithm 3 to compute (U k , ȳk ).The computational bottleneck of the procedure given in Algorithm 3 is the solution of the linear systems (4.13) and (4.14).Due to their Table 4.1 Jacobian's block times a vector: number of flops Algorithm 3 Nonlinear Gauss-Seidel algorithm input: y k−1 , U k−1 , µ k , η 2 from Algorithm 2 and ℓ max .

8:
end if 9: end for large dimensions we use a CG-like approach.The coefficient matrix in (4.13) takes the form: and it is positive semidefinite as Q k may be rank deficient.We can apply CG to (4.13) which is known to converge to the minimum norm solution if starting from the null approximation [21].Letting v ∈ R nr be the unitary eigenvector associated to the Since Q k does not depend on µ k and S k is expected to have r eigenvalues that go to zero as µ k , while the remaining ones stay bounded from above, we can conclude that the maximum eigenvalue of J T 11 J 11 + J T 21 J 21 does not depend on µ k while the smallest nonzero eigenvalue may go to zero at the same speed as µ 2 k .However in our numerical experimentation we verified that the term A T A acts as a regularization and the smallest nonzero eigenvalue of J T 11 J 11 + J T 21 J 21 remains bounded away from µ k also in the later iterations of the interior point algorithm.We will report on this later on, in the numerical results section (see Figure 6.1).
Let us now consider system (4.14).The coefficient matrix takes the form and it is positive definite.We repeat the reasoning applied earlier to J T 11 J 11 + J T 21 J 21 and conclude that Analogously we have Taking into account that r eigenvalues of X k do not depend on µ k while the remaining are equal to µ k , we conclude that the condition number of J T 22 J 22 increases as O(1/µ 2 k ).In the next subsection we will show how this matrix can be preconditioned.

Preconditioning J
Let us consider a preconditioner P k of the form with This choice is motivated by the fact that we discard the term 2µ k A(I n ⊗ U k U T k )A T from the expression of J T 22 J 22 .In fact, we use the approximation A similar idea is used in [33].An alternative choice involves matrix Z k of a smaller dimension This corresponds to introducing a further approximation We will analyze spectral properties of the matrix J T 22 J 22 preconditioned with P k defined in (4.19) with Z k given in (4.20).
Theorem 4.1.Let P k be given in (4.19) with Z k given in (4.20) and σ min (A) and σ max (A) denote the minimum and maximum singular values of A, respectively.The eigenvalues of the preconditioned matrix where ξ 1 and ξ 2 have the following forms: ) and 2  .
Proof.Note that Then, Let us denote with λ M and λ m the largest and the smallest eigenvalues of matrix , respectively.From (4.19) we deduce Then, using the Weyl inequality we obtain Moreover, Then, noting that .
Consequently, the eigenvalues of the preconditioned matrix belong to the interval (1 + ξ 1 , 1 + ξ 2 ), and the thesis follows.Note that from the result above, as µ k approaches zero, the minimum eigenvalue of the preconditioned matrix goes to one and the maximum remains bounded.
The application of P k to a vector d, needed at each CG iteration, can be performed through the solution of the (m + nr) × (m + nr) sparse augmented system: In order to recover the vector u = P −1 k d, we can solve the linear system and compute u as follows This process involves the inversion of AA T which can be done once at the beginning of the iterative process, and the solution of a linear system with matrix Note that E k has dimension n 2 × n 2 in case of choice (4.20) and dimension nr × nr in case of choice (4.21).Then, its inversion is impractical in case (4.20).On the other hand, using (4.21) we can approximately solve (4.23) using a CG-like solver.At this regard, observe that the entries of E k decay fast when far away from the main diagonal and E k can be preconditioned by its block-diagonal part, that is by where B is the operator that extracts from a matrix nr × nr its block diagonal part with n diagonal blocks of size r × r.

SDP reformulation of matrix completion problems.
We consider the problem of recovering a low-rank data matrix B ∈ R n×n from a sampling of its entries [12], that is the so called matrix completion problem.The problem can be stated as min rank( X) where Ω is the set of locations corresponding to the observed entries of B and the equality is meant element-wise, that is X s,t = B s,t , for all (s, t) ∈ Ω.Let m be the cardinality of Ω and r be the rank of B.
A popular convex relaxation of the problem [12] consists in finding the minimum nuclear norm of X that satisfies the linear constraints in (5.1), that is, solving the following heuristic optimization where the nuclear norm • * of X is defined as the sum of its singular values.
Candès and Recht proved in [12] that if Ω is sampled uniformly at random among all subset of cardinality m then with large probability, the unique solution to (5.2) is exactly B, provided that the number of samples obeys m ≥ C n5/4 r log n, for some positive numerical constant C. In other words, problem (5.2) is "formally equivalent" to problem (5.1).Let where X ∈ R n×n is the matrix to be recovered and W 1 , W 2 ∈ SR n×n .Then problem (5.2) can be stated as an SDP of the form (1.1) as follows where for each (s, t) ∈ Ω the matrix Θ st ∈ R n×n is defined element-wise for k, l = 1, . . ., n as see [28].We observe that primal variable X takes the form (5.3) with n = 2n, the symmetric matrix C in the objective of (1.1) is a scaled identity matrix of dimension n× n.The vector b ∈ R m is defined by the known elements of B and, for i = 1, . . ., m, each constraint matrix A i , corresponds to the known elements of B stored in b i .Matrices A i have a very special structure that yields nice properties in the packed matrix A. In fact, A ∈ R m×n 2 has 2m nonzero elements only and AA T = 1 2 I m and A(I n ) 2 = 0.
We now discuss the relationship between a rank r solution X of problem (5.2) and a rank r solution X of problem (5.4).
Proof.Let X = U ΣU T with U ∈ R 2n×r and Σ = R r×r be the singular value decomposition (SVD) of X.Let U be partitioned by that is X = U 1 ΣU T 2 has rank r.To prove the second part of the proposition, let X = U ΣV T with U, V ∈ R n×r and Σ = R r×r be the SVD factorization of X.We get the thesis by defining Assume that X has the form with U ∈ R n×r full column rank and µ ∈ R, then X has rank r.Proposition 5.3.If X is a rank r solution of (5.4), then X is a rank r solution of (5.2).Vice-versa, if X is a rank r solution of (5.2), then (5.4) admits a rank r solution.
Let X be a rank r solution of (5.2), t * = X * and U ΣV T , with U, V ∈ R n×r and Σ ∈ R r×r , be the SVD decomposition of X.Let us define X = W 1 X XT W 2 with W 1 = U ΣU T and W 2 = V ΣV T .Then X solves (5.4).In fact, X is positive semidefinite and 1  2 I • X = 1 2 (T race(W 1 ) + T race(W 2 )) = X * = t * .This implies that t * is the optimal value of (5.4).In fact, if we had Y such that Lemma 1] there would exist Ȳ such that Ȳ * < t * , that is Ȳ * < X * = t * .This is a contradiction as we assumed that t * is the optimal value of (5.2).
Remark.Assuming that a rank r solution to (5.2) exists, the above analysis justifies the application of our algorithm to search for a rank r solution of the SDP reformulation (5.4) of (5.2).We also observe that at each iteration our algorithm computes an approximation X k of the form Then, if at each iteration U k is full column rank, by Corollary 5.2, it follows that we generate a sequence { Xk } such that Xk has exactly rank r at each iteration k and it approaches a solution of (5.2).
Finally, let us observe that nnz(A) = m and m < n2 = n 2 /4.Then, by the analysis carried out in Subsection 4.1 each evaluation of the gradient of φ µ k amounts to O(n 2 r) flops and assuming to use a first-order method at each iteration to compute (U k , ȳk ) and to perform at most O(n) iterations of such method, each iteration of our method requires at most O(n 3 r) flops.

Numerical experiments on matrix completion problems.
We consider an application to matrix completion problems by solving (5.4) with our relaxed Interior Point algorithm for Low-Rank SDPs (IPLR), described in Algorithm 2. IPLR has been implemented using Matlab (R2018b) and all experiments have been carried out on Intel Core i5 CPU 1.3 GHz with 8 GB RAM.Parameters in Algorithm 2 have been chosen as follows: while the starting dual feasible approximation has been chosen as y 0 = 0, S 0 = 1 2 I n and U 0 is defined by the firs r columns of the identity matrix I n .
We perform two sets of experiments: the first aims at validating the proposed algorithms and is carried out on randomly generated problems; the second is an application of the new algorithms to a real data set.
6.1.Tests on random matrices.We generated n × n matrices of rank r by sampling two n × r factors B L and B R independently, each having independently and identically distributed Gaussian entries, and setting B = B L B R .The set of observed entries Ω is sampled uniformly at random among all sets of cardinality m.The matrix B is declared recovered if the (2,1) block X extracted from the solution X of (5.4), satisfies see [12].
Given r, we chose m by setting m = cr(2n − r), n = 600, 700, 800, 900, 1000.We used c = 0.01n + 4.These corresponding values of m are much lower than the theoretical bound provided by [12] and recalled in Section 5, but in our experiments they were enough to recover the sought matrix by IPLR.
In our experiments, the accuracy level in the matrix recovery in (6.1) is always achieved by setting ǫ = 10 −4 in Algorithm 2.
Let IPLR-GS denote the implementation of IPLR where the Gauss-Seidel strategy described in Algorithm 3 is used to perform Line 3 of Algorithm 2. We impose a maximum number of 5 ℓ-iterations and use the (possibly) preconditioned conjugate gradient method to solve the linear systems (4.13) and (4.14).We set a maximum of 100 CG iterations and the tolerance 10 −6 on the relative residual of the linear systems.System (4.13) is solved with unpreconditioned CG.Regarding (4.14), for the sake of comparison, we report statistics using unpreconditioned CG and CG employing the preconditioner defined by (4.19) and (4.21).In this latter case the action of the preconditioner has been implemented through the augmented system (4.22),following the procedure outlined at the end of Section 5.The linear system (4.23) has been solved by preconditioned CG, with preconditioner (4.24) allowing a maximum of 100 CG iterations and using a tolerance 10 −8 .In fact, the linear system (4.14)along the IPLR iterations becomes ill-conditioned and the application of the preconditioner needs to be performed with high accuracy.We will refer to the resulting method as IPLR-GS P.
As an alternative implementation to IPLR-GS, we considered the use of a firstorder approach to perform Line 3 of Algorithm 2. We implemented the Barzilai-Borwein method [3] with a non-monotone line-search following [15, Algorithm 1] and using parameter values as suggested therein.The Barzilai-Borwein method iterates until ∇φ µ k (U k , y k ) ≤ min(10 −3 , µ k ) or a maximum of 300 iterations is reached.We refer to the resulting implementation as IPLR-BB.
In the forthcoming tables we report: dimensions n and m of the resulting SDPs and target rank r of the matrix to be recovered; being X and S the computed solution, the final primal infeasibility A(X) − b , the complementarity gap XS − µI F , the error in the solution of the matrix completion problem E = X − B F / B F , the overall cpu time in seconds.
In Tables 6.1 and 6.2 we report statistics of IPLR-GS and IPLR-BB, respectively.We choose as a starting rank r the rank of the matrix B to be recovered.In the last column of Table 6.1 we report both the overall cpu time of IPLR-GS without preconditioner (cpu) and with preconditioner (cpu P) in the solution of (4.14).In bold the lowest computational time for each run.
As a first comment, we verified that Assumption 1 in Section 2 holds in our experiments.In fact, the method manages to preserve positive definiteness of the dual variable and α k < 1 is taken only in the early stage of the iterative process.Secondly, we observe that both IPLR-GS and IPLR-BB provide an approximation to the solution of the sought rank; in some runs the updating procedure increases the rank, but at the subsequent iteration the downdating strategy is activated and the procedure comes back to the starting rank r.Moreover, IPLR-GS is overall less expensive than IPLR-BB in terms of cpu time, in particular as n and m increase.In fact, the cost of the linear algebra in the IPLR-GS framework is contained as one/two inner Gauss-Seidel iterations are performed at each outer IPLR-GS iteration except for the very few initial iterations where up to five inner Gauss-Seidel iterations are needed.To give more details of the computational cost of both methods, in Table 6.3 we report some statistics of IPLR-GS and IPLR-BB for n = 900, r = 3 and 8.More precisely we report the average number of inner Gauss-Seidel iterations (avr GS) and the average number of unpreconditioned CG iterations in the solution of (4.13) (avr CG 1) and (4.14) (avr CG 2) for IPLR-GS and the average number of BB iterations for IPLR-BB (avr BB).We notice that the solution of SDP problems becomes more demanding as the rank increases, anyway both the number of BB iterations and the number of CG iterations are reasonable.To provide an insight into the linear algebra phase, in Figure 6.1 we plot the minimum nonzero eigenvalue and the maximum eigenvalue of the coefficient matrix of (4.13), i.e.

IPLR-GS
We remark that the matrix depends both Outer/Inner iterations on the outer iteration k and on the inner Gauss-Seidel iteration ℓ and we dropped the index ℓ to simplify the notation.Eigenvalues are plotted against the inner/outer iterations, for n = 100, r = 4 and IPLR-GS iterates until µ k < 10 −7 .In this run only one inner iteration is performed at each outer iteration except for the first outer iteration.We also plot in the left picture of Figure 6.2 the number of CG iterations versus inner/outer iterations.The figures show that the condition number of Q k and the overall behaviour of CG do not depend on µ k .Moreover, Table 6.3 shows that unpreconditioned CG is able to push the relative residual below 10 −6 in a low number of iterations even in the solution of larger problems and higher rank.These considerations motivate our choice of solving (4.13) without employing any preconditioner.
We now discuss the effectiveness of the preconditioner P k given in (4.19), with Zk given in (4.21), in the solution of (4.14).Considering n = 100, r = 4, in Figure 6.3 we plot the eigenvalue distribution (in percentage) of A(I ⊗ X 2 k )A T and P −1 k (A(I ⊗ X 2 k )A T ) at the first inner iteration of outer IPLR-GS iteration corresponding to µ k ≈ 1.9e−3.We again drop the index ℓ.We can observe that the condition number of the preconditioned matrix is about 1.3e5, and it is significantly smaller than the condition number of the original matrix (about 3.3e10).The preconditioner succeeded both in pushing the smallest eigenvalue away from zero and in reducing the largest eigenvalue.However, CG converges in a reasonable number of iterations even in the unpreconditioned case, despite the large condition number.In particular, we can observe in the right picture of Figure 6.2 that preconditioned CG takes less than five iterations in the last stages of IPLR-GS and that the most effort is made in the initial stage of the IPLR-GS method; in this phase the preconditioner is really effective in reducing the number of CG iterations.These considerations remain true even for larger values of n and r as it is shown in is Table 6.3.
Focusing on the computational cost of the preconditioner's application, we can observe from the cpu times reported in Table 6.1, that for r = 3, 4, 5 the employment of the preconditioner produces a great benefit, with savings that vary from 20% to 50%.Then, the overhead associated to the construction and application of the preconditioner is more than compensated by the gains in the number of CG iterations.The cost of application of the preconditioner increases with r as the dimension of the diagonal blocks of M k in (4.24) increases with r.Then, for small value of n and r = 6, 7, 8 unpreconditioned CG is preferable, while for larger value of n the preconditioner is effective in reducing the overall computational time for r ≤ 7.This behaviour is summarized in Figure 6.4 where we plot the ratio cpu P/cpu with respect to dimension n and rank (from 3 to 8).
In the approach proposed in this paper the primal feasibility is gradually reached, hence it is also possible to handle data B Ω corrupted by noise.To test how the method behaves in such situations we set B(s,t) = B (s,t) + ηRD (s,t) for any (s, t) ∈ Ω, where RD (s,t) is a random scalar drawn from the standard normal distribution, generated by the Matlab function randn; η > 0 is the level of noise.Then, we solved problem (5.4) using the corrupted data B(s,t) to form the vector b.Note that, in this case In order to take into account the presence of noise we set ǫ = max(10 −4 , 10 −1 η) in Algorithm 2.
Results of these runs are collected in Table 6.4 where we considered η = 10 −1 and started with the target rank r.In the table we report X − B F /n, that is the root-mean squared error per entry.Note that the root-mean error per entry in data B Ω is of the order of the noise level 10 −1 , as well as A(B) − b 2 / √ m.Then, we claim to recover the matrix with acceptable accuracy, corresponding to an average error smaller than the level of noise.Rank updating.We now test the effectiveness of the rank updating/downdating strategy described in Algorithm 2. To this purpose, we run IPLR-GS P starting from r = 1 and report the results in Table 6.5 for n = 600, 800, 1000.In all runs, the target rank has been correctly identified by the updating strategy and the matrix B is well-recovered.Runs in italic have been obtained allowing 10 inner Gauss-Seidel iterations.In fact, 5 inner Gauss-Seidel iterations were not enough to sufficiently reduce the residual in (2.4) and the procedure did not terminate with the correct rank.Comparing the values of the cpu time in Tables 6.1 and 6.5 we observe that the cities in the US and Canada computed from latitude/longitude data.

IPLR-GS
We sampled the 30% of the matrix G of geodesic distances and computed a lowrank approximation X by IPLR-BB inhibiting rank updating/downdating and using ǫ = 10 −4 .We compared the obtained approximation with the best rank-r approximation Xr , computed by truncated SVD, that requires the knowledge of the full matrix G.We considered some small values for the rank (r = 3, 4, 5) and reported in Table 6.7 the errors E = G − X F / G F and E r = G − Xr F / G F .We remark that the matrix G is not nearly-low-rank, and the method correctly detects that there does not exist a feasible rank r matrix as it is not able to decrease the primal infeasibility below 1e0.On the other hand the error in the provided approximation, obtained using only the 23% of the entries, is the same as that of the best rank-r approximation Xr .Note that computing the 5-rank approximation is more demanding.In fact the method requires on average: 3.4 Gauss-Seidel iterations, 37 unpreconditioned CG iterations for computing ∆U and 18 preconditioned CG iterations for computing ∆y.In contrast, the 3-rank approximation requires on average: 3.8 Gauss-Seidel iterations, 18 unpreconditioned CG iterations for computing ∆U and 10 preconditioned CG iterations for computing ∆y.As a final comment, we observe that IPLR-GS fails when r = 5 since unpreconditioned CG struggles with the solution of (4.14).The computed direction ∆y is not accurate enough and the method fails to maintain S positive definite within the maximum number of allowed backtracks.Applying the preconditioner cures the problem because more accurate directions become available.

Conclusions.
We have presented a new framework for an interior point method for low rank semidefinite programming.The method relaxes the rigid IPM structure and replaces the general matrix X with the special form (2.3) which by construction enforces a convergence to a low rank solution as µ goes to zero.Therefore effectively instead of requiring a general n×n object, the proposed method works with an n × r matrix U , which delivers significant storage and cpu time savings.It also handles well problems with noisy data and allows to adaptively correct the (unknown) rank.We have analyzed the convergence of the method and have proposed an efficient implementation which accommodates both the first-and the second-order methods to compute search directions.Tests have been performed using Barzilai-Borwein and Gauss-Seidel methods, respectively.Extensive computational evidence has been provided to demonstrate the efficiency of the proposed method and its ability to handle large scale matrix completion problems.Regarding the choice of first-or second-order methods to compute search directions, in our experience both approaches performed well.The Barzilai-Borwein approach is slower than the Gauss-Seidel, but it is simpler
Fig.6.4.The ratio cpu P/cpu as a funcion of dimension n and of the rank (data extracted from Table6.1).

Table 6 .
P 7 IPLR-GS P for low rank approximation of the City Distance matrix.