A Rank-Corrected Procedure for Matrix Completion with Fixed Basis Coefficients

Introduction

The low-rank matrix completion is to recover an unknown low-rank matrix from the under-sampled observations with or without noises. This problem is of considerable interest in many application areas, from machine learning to quantum state tomography. A basic idea to address a low-rank matrix completion problem is to minimize the rank of a matrix subject to certain constraints from observations. Since the direct minimization of rank function is generally NP-hard, a widely-used convex relaxation approach is to replace the rank function with the nuclear norm — the convex envelope of the rank function over a unit ball of the spectral norm .

The nuclear norm technique has been observed to provide a low-rank solution in practice for a long time (see, e.g., ). The first remarkable theoretical characterization for the minimum rank solution via the nuclear norm minimization was given by Recht, Fazel and Parrilo , with the help of the concept of Restricted Isometric Property (RIP). Recognizing that the matrix completion problem does not obey the RIP, Candès and Recht introduced the concept of incoherence property and proved that most low-rank matrices can be exactly recovered from a surprisingly small number of noiseless observations of randomly sampled entries via the nuclear norm minimization. The bound of the number of sampled entries was later improved to be near-optimal by Candès and Tao through a counting argument. Such a bound was also obtained by Keshavan et al. for their proposed OptSpace algorithm. Later, Gross sharpened the bound by employing a novel technique from quantum information theory developed in , in which noiseless observations were extended from entries to coefficients relative to an arbitrary basis. This technique was also adapted by Recht , leading to a short and intelligible analysis. Besides the above results for the noiseless case, matrix completion with noise was first addressed by Candès and Plan . More recently, nuclear norm penalized estimators for matrix completion with noise have been well studied by Koltchinskii, Lounici and Tsybakov , Negahban and Wainwright , and Klopp under different settings. Besides the nuclear norm, estimators with other penalties for matrix completion have also been considered in terms of recoverability in the literature, e.g., .

The nuclear norm technique has been demonstrated to be a successful approach to encourage a low-rank solution for matrix completion. However, its efficiency may be challenged in some circumstances. For example, Salakhutdinov and Srebro showed that when certain rows and/or columns are sampled with high probability, the nuclear norm minimization may fail in the sense that the number of observations required for recovery is much more than the setting of most matrix completion problems. It means that the efficiency of the nuclear norm techniques could be highly weakened under a general sampling scheme. Negahban and Wainwright also pointed out the impact of such heavy sampling schemes on the recovery error bound. As a remedy for this, a weighted nuclear norm (trace norm), based on row- and column-marginals of the sampling distribution, was suggested in if the prior information on sampling distribution is available. Moreover, the conditions characterized by Bach for rank consistency of the nuclear norm penalized least squares estimator may not be satisfied, especially when certain constraints are involved.

A concrete example of interest is to recover a density matrix of a quantum system from Pauli measurements in quantum state tomography (see, e.g., ). A density matrix is a Hermitian positive semidefinite matrix of trace one. Clearly, if the constraints of positive semidefiniteness and trace one are simultaneously imposed on the nuclear norm minimization, the nuclear norm completely fails in promoting a low-rank solution. Thus, one of the two constraints has to be abandoned in the nuclear norm minimization and then be restored in the post-processing stage. In fact, this idea has been much explored in and the numerical results there indicated its relative efficiency though it still has much room for improvement.

All the above examples motivate us to ask whether it is possible to go beyond the nuclear norm approach for practical use to seek for better performance in low-rank matrix completion. In this paper, we provide a positive answer to this question with both theoretical and empirical supports. We first establish a unified low-rank matrix completion model, which allows for the imposition of fixed basis coefficients so that the correlation and the density matrix completion are included as special cases. It means that in our setting, for any given basis of the matrix space, a few basis coefficients of the true matrix are assumed to be fixed due to a certain structure or some prior information, and the rest are allowed to be observed with noises under a general sampling scheme. To pursue a low-rank solution with a high recovery accuracy, we propose a rank-correction step to generate a new estimator. The rank-correction step solves a penalized least squares problem with its penalization being the nuclear norm minus a linear rank-correction term constructed on a reasonable initial estimator. A satisfactory choice of the initial estimator could be the nuclear norm penalized least squares estimator or one of its analogies. The resulting convex matrix optimization problem can be solved by the efficient algorithms recently developed in even for large-scale cases.

The idea of using a two-stage or even multi-stage procedure is not brand new for dealing with sparse recovery in the statistical and machine learning literature. The $l_{1}$-norm penalized least squares method, also known as the Lasso , is very attractive and popular for variable selection in statistics, thanks to the invention of the fast and efficient LARS algorithm . On the other hand, the $l_{1}$-norm penalty has long been known by statisticians to yield biased estimators and cannot achieve the best estimation performance . The issue of bias can be overcome by nonconvex penalization methods, see, e.g., . A multi-stage procedure naturally occurs if the nonconvex problem obtained is solved by an iterative algorithm . In particular, once a good initial estimator is used, a two-stage estimator is enough to achieve the desired asymptotic efficiency, e.g., the adaptive Lasso proposed by Zou . There are also a number of important works along this line on variable selection, including , to name only a few. For a broad overview, the interested readers are referred to the recent survey papers . It is natural to extend the ideas from the vector case to the matrix case. Fazel, Hindi and Boyd first proposed the reweighted trace minimization for minimizing the rank of a positive semidefinite matrix. In , Bach made an important step in extending the adaptive Lasso of Zou to the matrix case for rank consistency. However, it is not clear how to apply Bach’s idea to our matrix completion model with fixed basis coefficients since the required rate of convergence of the initial estimator for achieving asymptotic properties is no longer valid, as far as we can see. More critically, there are numerical difficulties in efficiently solving the resulting optimization problems. Numerical difficulties also occur in the reweighted nuclear norm approach proposed by Mohan and Fazel as an extension of for rectangular matrices. Iterative reweighted least squares minimization is an alternative extension of independently proposed by Mohan and Fazel and Fornasier, Rauhut and Ward , taking advantage of the property that the rank of a matrix is equal to the rank of the product of this matrix and its transpose. However, the resulting smoothness of inner-iteration subproblems is weak in encouraging a low-rank solution so much more iterations are needed in general and thus the computational cost is high especially when hard constraints such as fixed basis coefficients are involved.

The rank-correction step to be proposed in this paper is for overcoming the above difficulties. This approach is inspired by the majorized penalty method proposed by Gao and Sun for solving structured matrix optimization problems with a low-rank constraint. For our proposed rank-correction step, we establish a non-asymptotic recovery error bound in Frobenius norm, following a similar argument adopted by Klopp in . We also discuss the impact of adding the rank-correction term on recovery error. More importantly, we provide an affirmative guarantee that under mild condition the rank-correction step highly improves the recoverability, compared with the nuclear norm penalized least squares estimator. As the estimator is expected to be of low-rank, we also study the asymptotic property — rank consistency in the sense of Bach , under the setting that the matrix size is assumed to be fixed. This setting may not be ideal for analyzing asymptotic properties for matrix completion, but it does allow us to take the crucial first step to gain insights into the limitation of the nuclear norm penalization. Among others, the concept of constraint nondegeneracy for conic optimization problem plays a key role in our analysis. Interestingly, our results of recovery error bound and rank consistency suggest a consistent criterion for constructing a suitable rank-correction function. In particular, for the correlation and the density matrix completion problems, we prove that rank consistency automatically holds for a broad selection of rank-correction functions. For most cases, a single rank-correction step is sufficient for a substantial improvement, unless the sample ratio is rather low so that the rank-correction step may be iteratively used for two or three times to achieve the limit of improvement. Owing to this property, the advantage of our proposed method is more apparent in practical computations especially when fixed basis coefficients are involved. Finally, we remark that our results can also be used to provide a theoretical foundation in the statistical setting for the majorized penalty method of Gao and Sun and Gao for structured low-rank matrix optimization problems.

This paper is organized as follows. In Section 2, we introduce the observation model of matrix completion with fixed basis coefficients and formulate the rank-correction step. In Section 3, we establish a non-asymptotic recovery error bound for the estimator generated from the rank-correction step and provide a quantification of the improvement in recoverability. Section 4 provides necessary and sufficient conditions for rank consistency. Section 5 is devoted to the construction of the rank-correction function. In Section 6, we report numerical results to validate the efficiency of our proposed rank-corrected procedure. We conclude this paper in Section 7. All relevant material and all proofs of theorems are left in the appendices.

Notation. Here we provide a brief summary of the notation used in this paper.

The notation T denotes the transpose for the real case and the conjugate transpose for the complex case. The notation ∗ means the adjoint of a linear operator.

For any given vector $x$, $\text{Diag}(x)$ denotes a rectangular diagonal matrix of suitable size with the $i$-th diagonal entry being $x_{i}$.

The notations $\stackrel{{\scriptstyle a.s.}}{{\rightarrow}}$, $\stackrel{{\scriptstyle p}}{{\rightarrow}}$ and $\stackrel{{\scriptstyle d}}{{\rightarrow}}$ mean almost sure convergence, convergence in probability and convergence in distribution, respectively. We write $x_{m}=O_{p}(1)$ if $x_{m}$ is bounded in probability.

For any set $K$, let $\delta_{K}(x)$ denote the indicator function of $K$, i.e., $\delta_{K}(x)=0$ if $x\in K$, and $\delta_{K}(x)=+\infty$ otherwise. Let $I_{n}$ denote the $n\times n$ identity matrix.

Problem formulation

In this section, we formulate the model of the matrix completion problem with fixed basis coefficients, and then propose an adaptive nuclear semi-norm penalized least squares estimator for solving this class of problems.

When a few basis coefficients are fixed, one only needs to observe the rest for recovering the unknown matrix $\overline{X}$. Assume that we are given a collection of $m$ noisy observations of the basis coefficients relative to $\{\Theta_{k}:k\in\beta\}$ in the following form

The indices $\omega_{1},\ldots,\omega_{m}$ are i.i.d. copies of a random variable $\omega$ that has a probability distribution $\Pi$ over $\{1,\ldots,d\}$ defined by

Note that each $\Theta_{k},k\in\beta$ is assumed to be sampled with a positive probability in this sampling scheme. In particular, when the sampling probability of all $k\in\beta$ are equal, i.e., $p_{k}=1/d_{2}\ \forall\,k\in\beta$, we say that the observations are sampled uniformly at random.

Then, the observation model (1) can be expressed in the following vector form

Next, we present some examples of low-rank matrix completion problems in the above settings.

Here, $\sqrt{-1}$ represents the imaginary unit. Of course, one may fix some off-diagonal entries in specific applications.

Density matrix completion. A density matrix of dimension $n=2^{l}$ for some positive integer $l$ is an $n\times n$ Hermitian positive semidefinite matrix with trace one. In quantum state tomography, one aims to recover a density matrix from Pauli measurements (observations of the coefficients relative to the Pauli basis) , given by

where “$\otimes$” means the Kronecker product of two matrices and

2 The rank-correction step

In many situations, the nuclear norm penalization performs well for matrix recovery, but its efficiency may be challenged if the observations are sampled at random obeying a general distribution such as the one considered in . The setting of fixed basis coefficients in our matrix completion model can also be regarded to be under an extreme sampling scheme. In particular, for the correlation and density matrix completion, the nuclear norm completely loses its efficiency since it reduces to a constant in these two cases. In order to overcome the shortcomings of the nuclear norm penalization, we propose a rank-correction step to generate an estimator in pursuit of a better recovery performance.

Recall that $\overline{X}$ is the unknown true matrix of rank $r$. Given an initial estimator $\widetilde{X}_{m}$ of $\overline{X}$, say, the nuclear norm penalized least squares estimator or one of its analogies, our proposed rank-correction step is to solve the convex optimization problem

Hereafter, we call $F$ the rank-correction function and $\langle F(\widetilde{X}_{m}),X\rangle$ the rank-correction term. Note that, when $F\equiv 0$, the rank-correction step (3) reduces to the nuclear norm penalized least squares estimator, which equally penalizes singular values to promote a low-rank solution for matrix completion. Certainly, for this purpose, penalizing more on small singular values or even directly penalizing the rank function could serve better, but only theoretically rather than practically, due to the lack of convexity. Also note that an initial estimation, if deviates not too much from the true matrix, could contain some information of the singular values and/or the rank of the true matrix to a certain extent. Therefore, provided such an initial estimator is available, it is achievable to construct a rank-correction term with a suitable $F$ to substantially offset the penalization of large singular values from the nuclear norm penalty. Consequently, we can expect the rank-correction step (3) to have a better low-rank promoting ability and outperform the nuclear norm penalized least squares estimator.

The key issue is then how to construct a favored rank-correction function $F$. In the next two sections, we provide theoretical supports to our proposed rank-correction step, from which some important guidelines on the construction of $F$ can be captured. In particular, if one chooses the nuclear norm penalized least squares estimator to be the initial estimator $\widetilde{X}_{m}$, and also suitably chooses the spectral operator $F$ so that $\|X\|_{*}-\langle F(\widetilde{X}_{m}),X\rangle$ is a semi-norm, called nuclear semi-norm, then the estimator $\widehat{X}_{m}$ generated from this two-stage procedure is called the adaptive nuclear semi-norm penalized least squares estimator associated with $F$.

3 Relation with the majorized penalty approach

The rank-correction step above is inspired by the majorized penalty approach proposed by Gao and Sun for solving the rank constrained matrix optimization problem:

where $\|X\|_{(r)}:=\sigma_{1}(X)+\cdots+\sigma_{r}(X)$ denotes the Ky Fan $r$-norm. The central idea of the majorized penalty approach is to solve the following penalized version of (4):

where $\rho>0$ is the penalty parameter. With the current iterate $X^{k}$, the majorized penalty approach yields the next iterate $X^{k+1}$ by solving the convex optimization problem

where $G^{k}$ is a subgradient of the convex function $\|X\|_{(r)}$ at $X^{k}$, and $\widehat{h}^{k}$ is a convex majorization function of $h$ at $X^{k}$. By comparing with (3), one may notice that our proposed rank-correction step is close to a single step of the majorized penalty approach.

Note that the rank constrained least squares problem is of great consideration in matrix completion especially when the rank information is known. However, different from the noiseless case, for matrix completion with noise, the solution to the rank constrained least squares problem (assuming the uniqueness) is in general not the true matrix though quite close to it. Indeed, there may exist many candidate matrices surrounding the true matrix and having its rank. The rank constrained least squares solution is only one of them. It deviates the least from the noisy observations rather than the true matrix. Naturally, it is conceivable that some candidate matrices may deviate a bit more from the noisy observations but less from the true matrix. So, for the purpose of matrix completion, there is no need to aim precisely at the rank constrained least squares solution and find this solution accurately. An approach roughly towards it such as our proposed rank-correction step (3) is good enough to bring similar good recovery performance.

Error bounds

In this section, we aim to derive a recovery error bound in Frobenius norm for the estimator generated from the rank-correction step (3) and discuss the impact of the rank-correction term on the resulting bound. The analysis mainly follows Klopp’s arguments in , which is also in line with those used by Negahban and Wainwright .

We start the analysis by defining a quantity, which plays a key role in the subsequent analysis, as

A basic relation between the true matrix $\overline{X}$ and its estimate $\widehat{X}_{m}$ can be obtained by using the optimality of $\widehat{X}_{m}$ to the problem (3) as follows.

For any $\kappa>1$, if \rho_{m}\geq\kappa\nu\Big{\|}\frac{1}{m}\mathcal{R}_{\Omega}^{*}(\xi)\Big{\|}, then the following inequality holds:

We emphasize that $\kappa$ is not restricted to be a constant in Theorem 1 but could be set to depend on the size of matrix. This realization is important as can be seen in the sequel. According to Theorem 1, the choice of the penalty parameter $\rho_{m}$ depends on the observation noises $\xi_{i}$ and the sampling operator $\mathcal{R}_{\Omega}$. Therefore, we make the following assumption on the noises $\xi_{i}$ as follows:

The i.i.d. noise variables $\xi_{i}$ are sub-exponential, i.e., there exist positive constants $c_{1}$, $c_{2}$ and $c_{3}$ such that for all $t>0$, ${\rm Pr}(|\xi_{i}|\geq t)\leq c_{1}\exp(-c_{2}t^{c_{3}}).$

Moreover, based on Assumption 1, we further define quantities $\mu_{1}$ and $\mu_{2}$ that control the sampling probability for observations as

Theorem 1 reveals the key to deriving a recovery error bound in Frobenius norm, that is, to establish the relation between $\frac{1}{m}\|\mathcal{R}_{\Omega}(\widehat{X}_{m}-\overline{X})\|_{2}^{2}$ and $\|\widehat{X}_{m}-\overline{X}\|_{F}^{2}$. This can be achieved by looking into some RIP-like property of the sampling operator $\mathcal{R}_{\Omega}$, as done previously in . Following this idea, we obtain an explicit recovery error bound as follows:

Under Assumptions 1 and 2, there exist some positive absolute constants $c_{0},c_{1},c_{2},c_{3}$ and some positive constants $C_{0},C_{1}$ (only depending on the $\psi_{1}$ Orlicz norm of $\xi_{k}$) such that when $m\geq c_{3}\sqrt{d_{2}}\log^{3}(n_{1}+n_{2})/\mu_{2}$, for any $\kappa>1$, if $\rho_{m}$ is chosen as

then with probability at least $1-c_{1}(n_{1}+n_{2})^{-c_{2}}$,

Theorem 2 shows that for any rank-correction function $F$, controlling the recovery error only needs the samples size $m$ to be of roughly the degree of freedom of a rank $r$ matrix up to a logarithmic factor in the matrix size. Besides the information on the order of magnitude, Theorem 2 also provides us more details on the constant part in the recovery error bound, which also plays an important role in practice. The impact of different choices of rank-correction functions on recovery error is fully embodied with the value of $a_{m}$. Note that the smaller $a_{m}$ is, the smaller the error bound (10) is for a fixed $\kappa$, and thus the smaller value this error bound can achieve for the best $\kappa$ (as well as the best $\rho_{m}$). Therefore, we aim to establish an explicit relationship between $a_{m}$ and $F$ in the next theorem.

where $\tau>0$. As can be seen, the larger the matrix size $n$ is, the easier $a_{m}$ becomes less than $1$ or even close to . If the rank of the true matrix is unknown, one could construct the rank-correction function $F$ on account of the tradeoff between optimality and robustness, to be discussed in Section 5. An experimental example of the relationship between $a_{m}$ and $F$ can be found in Table 1.

Next, we demonstrate the power of the rank-correction term with more details. It is interesting to notice that the value of $\kappa$ (as well as $\rho_{m}$) has a substantial impact on the recovery error bound (10). The part related to the magnitude of noise $\nu$ increases as $\kappa$ increases, while the part related to the upper bound $b$ of entries slightly decreases to its limit as $\kappa$ increases. Therefore, our first target is to find the smallest error bound in terms of (10) among all possible $\kappa>1$. It is possible to work on the error bound (10) directly for its minimum in $\kappa$ but the subsequent analysis is much more tedious. For simplicity of illustration, instead, we perform our analysis on a slightly relaxed version instead as

Direct calculation shows that over $\kappa>1$, $\eta_{m}$ attains its minimum

It is worthwhile to note that \overline{\kappa}=O\big{(}1/\sqrt{a_{m}}\,\big{)} when $a_{m}\ll 1$, meaning that the optimal choice of $\kappa$ is inversely proportional to $\sqrt{a_{m}}$ rather than a simple constant. (This observation is important for achieving the rank consistency in Section 4.) In other words, for achieving the best possible recovery error, the penalty parameter $\rho_{m}$ chosen for the rank-correction step (3) with $a_{m}<1$ should be larger than that for the nuclear norm penalized least squares estimator. In addition, consider two extreme cases with $a_{m}=1$ and $a_{m}=0$ respectively:

By direct calculations, we obtain $\overline{\eta}^{0}/\overline{\eta}^{1}\in(0.356,0.586)$, where the lower bound is attained when $c_{0}\nu=b$ and the upper bound is approached when $c_{0}\nu/b\rightarrow 0$ or $c_{0}\nu/b\rightarrow\infty$. This finding motivates us to wonder whether the recovery error can be reduced by around half in practice. This inference is further validated by numerical experiments in Section 6.

Rank consistency

In this section we consider the asymptotic behavior of the estimator generated from the rank-correction step (3) in term of its rank. We expect that the resulting $\widehat{X}_{m}$ has the same rank as the true matrix $\overline{X}$. Theorem 2 only reveals a flavored parameter $\rho_{m}$ in terms of the optimal order but rather its exact value. In practice, for a chosen parameter $\rho_{m}$, there is hardly any clue to know the recovery performance of the resulting solution since the true matrix is unknown. However, if the rank property holds as expected, the observable rank information may be used to infer the recovery quality of the resulting solution of a parameter and thus help in parameter searching. Numerical experiments in Section 6 demonstrate the practicability of this idea.

For the purpose above, we study the rank consistency in the sense of Bach under the setting that the matrix size is fixed. An estimator $X_{m}$ of the true matrix $\overline{X}$ is said to be rank consistent if

Throughout this section, we make the following assumptions:

The spectral operator $F$ is continuous at $\overline{X}$.

The initial estimator $\widetilde{X}_{m}$ satisfies $\widetilde{X}_{m}\stackrel{{\scriptstyle p}}{{\rightarrow}}\overline{X}$ as $m\rightarrow\infty$.

Epi-convergence in distribution gives us an elegant way in analyzing the asymptotic behavior of optimal solutions of a sequence of constrained optimization problems. Based on this technique, we obtain the following result.

If $\rho_{m}\rightarrow 0$, then $\widehat{X}_{m}\stackrel{{\scriptstyle p}}{{\rightarrow}}\overline{X}$ as $m\rightarrow\infty$.

For notational simplicity, we divide the index set $\beta$ into three subsets as

By extending the arguments of Bach for the nuclear norm penalized least squares estimator from the unconstrained case to the constrained case, we obtain the following results.

If $\rho_{m}\rightarrow 0$ and $\sqrt{m}\rho_{m}\rightarrow\infty$, then for the rank consistency of $\widehat{X}_{m}$,

For the positive semidefinite case, the nuclear norm $\|X\|_{*}$ in (3) simply reduces to the trace $\langle I_{n},X\rangle$. We assume that the Slater condition holds.

Under Assumption 5, if $\rho_{m}\rightarrow 0$ and $\sqrt{m}\rho_{m}\rightarrow\infty$, then for the rank consistency of $\widehat{X}_{m}$,

Clearly, the invertibility of $\mathcal{B}_{2}$ is equivalent to the uniqueness of the solution to the linear systems (13) and (14). The following result provides a link between the constraint nondegeneracy and the positive definiteness of $\mathcal{B}_{2}$.

Combining Theorems 5, 6 and 7 together with (18) and (19), we immediately have the following result of rank consistency.

Suppose that $\rho_{m}\rightarrow 0$ and $\sqrt{m}\rho_{m}\rightarrow\infty$. If

then the estimator $\widehat{X}_{m}$ generated from the rank-correction step (3) is rank consistent.

The covariance matrix completion with partial positive diagonal entries fixed.

Due to the positive semidefinite structure, the magnitudes of off-diagonal entries are fully controlled by the magnitudes of diagonal entries. Therefore, we remove all the bounded constraints corresponding to off-diagonal entries from the rank-correction step (3) as they are redundant. Thus, the constraints are reduced to

where $(\pi,\pi^{c})$ is a partition of the index set $\{1,\ldots,n\}$. This class of problems includes the correlation matrix completion as a special case, in which all diagonal entries are fixed to be ones.

The density matrix completion with its trace fixed to be one.

Due to the positive semidefinite structure, all the coefficients of Pauli basis are controlled because of the trace one constraint. Therefore, we remove all the bounded constraints from the rank-correction step (3) as they are redundant. Thus, in this case the constraints are reduced to

Interestingly, for the matrix completion problems of Classes I and II, the constraint nondegeneracy automatically holds at $\overline{X}$. More importantly, if observations are sampled uniformly at random, the rank consistency can be guaranteed for a broad class of rank-correction functions $F$.

then the estimator $\widehat{X}_{m}$ generated from the rank-correction step (3) is rank consistent.

Construction of the rank-correction function

Next, we proceed with the construction of the rank-correction function $F$ for the rectangular case. For the positive semidefinite case, one only needs to replace the singular value decomposition with the eigenvalue decomposition and conduct exactly the same analysis.

If the rank of the true matrix $\overline{X}$ is known, it is clear that the best choice of $F$ is

2 The rank is unknown

If the rank of the true matrix $\overline{X}$ is unknown, we intend to construct a spectral operator $F$ to imitate the case when the rank is known. Here, we propose $F$ to be a spectral operator

Let $F$ be a spectral operator defined by (24), (25) and (26).

If \frac{\|\widetilde{X}_{m}-\overline{X}\|_{F}}{\sigma_{r}(\overline{X})}<\frac{1}{\sqrt{2}}\big{(}1-e^{-\sqrt{2r}}\big{)}, then for any $\varepsilon$ satisfying $\frac{\sigma_{r+1}(\widetilde{X}_{m})}{\sigma_{1}(\widetilde{X}_{m})}<\varepsilon<\frac{\sigma_{r}(\widetilde{X}_{m})}{\sigma_{1}(\widetilde{X}_{m})}$, there exists some $\overline{\tau}_{1}>0$ such that $a_{m}<1$ for any $F$ with $\tau\geq\overline{\tau}_{1}$.

Suppose that the constraint nondegeneracy holds at $\overline{X}$ to the problem (3). If $\rho_{m}\rightarrow 0$ and $\sqrt{m}\rho_{m}\rightarrow\infty$, then for any $\varepsilon$ satisfying $0<\varepsilon<\frac{\sigma_{r}(\overline{X})}{\sigma_{1}(\overline{X})}$, there exists some $\overline{\tau}_{2}>0$ such that the rank consistency of $\widehat{X}_{m}$ holds for any $F$ with $\tau\geq\overline{\tau}_{2}$.

The proof of Corollary 10 is straightforward so we omit it. Corollary 10 suggests an ideal choice of $\varepsilon$ for the recovery error reduction, i.e., $\varepsilon\in\left(\frac{\sigma_{r+1}(\widetilde{X}_{m})}{\sigma_{1}(\widetilde{X}_{m})},\frac{\sigma_{r}(\widetilde{X}_{m})}{\sigma_{1}(\widetilde{X}_{m})}\right)$, provided that $\widetilde{X}_{m}$ does not deviate too much from $\overline{X}_{m}$, and also an ideal choice of $\varepsilon$ for rank consistency, i.e., $\varepsilon\in\left(0,\frac{\sigma_{r}(\overline{X}_{m})}{\sigma_{1}(\overline{X}_{m})}\right)$. Note that these two intervals may not overlap each other, implying the theoretical possibility that the recovery error reduction and the rank consistency may not be achieved simultaneously if the initial estimator $\widetilde{X}_{m}$ is not close to $\overline{X}_{m}$.

The interval of $\varepsilon$ for the recovery error reduction is disclosed if the true rank is accessible. Therefore, this ideal interval is an important insight that can be used to guide the choice of $\varepsilon$ in practice since the initial $\widetilde{X}_{m}$ should contain some information of the true rank in general. Indeed, the value of $\varepsilon$ can be regarded as a divide of confidence on whether $\sigma_{i}(\widetilde{X}_{m})$ is believed to come from a nonzero singular values of $\overline{X}$ with perturbation — positive confidence if $\sigma_{i}(\widetilde{X}_{m})>\varepsilon\sigma_{1}(\widetilde{X}_{m})$ and negative confidence if $\sigma_{i}(\widetilde{X}_{m})<\varepsilon\sigma_{1}(\widetilde{X}_{m})$. Next we look for a suitable $\tau$. It is observed from Figure 1 that the parameter $\tau>0$ mainly controls the shape of $\phi$ over $t\in$. The function $\phi$ is concave if $0<\tau\leq 1$ and $S$-shaped with a single inflection point at \varepsilon\big{(}\frac{\tau-1}{\tau+1}\big{)}^{1/\tau} if $\tau>1$. It should be good to choose an $S$-shaped function $\phi$. But one also needs to take account of the steepness of $\phi$, which increases when $\tau$ increases. In particular for any $\varepsilon$ satisfying $0<\varepsilon<1$, $\phi$ approaches to the step function taking the value if $0\leq t<\varepsilon$ and the value $1$ if $\varepsilon<t\leq 1$ as $\tau\rightarrow\infty$. Since the rank of $\overline{X}$ is unknown and the singular values of $\widetilde{X}_{m}$ are unpredictable, choosing a large $\tau$ could be risky. Therefore, one needs to choose $\tau$ with certain conservation, sacrificing certain recovery quality in exchange for robustness strategically. Here, we provide a recommendation of the choices $\varepsilon\approx 0.05$ (or within $0.01\sim 0.1$) and $\tau=2$ (or within $1\sim 3$) for most cases, particularly when the initial estimator is generated from the nuclear norm penalized least squares problem. These choices have performed very stably for plenty of problems, as validated in Section 6.

We also remark that for the positive semidefinite case, the rank-correction function defined by (24), (25) and (26) is related to the reweighted trace norm for the matrix rank minimization proposed by Fazel et al. . The reweighted trace norm in for the positive semidefinite case is $\langle(X^{k}+\varepsilon I_{n})^{-1},X\rangle$, which arises from the derivative of the surrogate function $\log\det(X+\varepsilon I_{n})$ of the rank at an iterate $X^{k}$, where $\varepsilon$ is a small positive constant. Meanwhile, in our proposed rank-correction step, if we choose $\tau=1$, then $I_{n}-\frac{1}{1+\varepsilon}F(\widetilde{X}_{m})=\varepsilon^{\prime}(\widetilde{X}_{m}+\varepsilon^{\prime}I_{n})^{-1}$ with $\varepsilon^{\prime}=\varepsilon\|\widetilde{X}_{m}\|$. Superficially, similarity occurs; however, it is notable that $\varepsilon^{\prime}$ depends on $\widetilde{X}_{m}$, which is different from the constant $\varepsilon$ in . More broadly speaking, the rank-correction function $F$ defined by (24), (25) and (26) is not a gradient of any real-valued function. This distinguishes our proposed rank-correction step from the reweighted trace norm minimization in even for the positive semidefinite case.

Numerical experiments

In this section, we validate the power of our proposed rank-correction step on the recovery by applying it to different matrix completion problems. We adopted the proximal alternating direction method of multipliers (proximal ADMM) to solve the optimization problem (3). For more details of the proximal ADMM, the readers may refer to Appendix B of . For convenience, in the sequel, the NNPLS estimator and the RCS estimator, respectively, stand for the estimators from the nuclear norm penalized least squares problem (i.e., $F\equiv 0$) and the rank-correction step (3) with $F$ specified in Section 5. Given an estimator $X_{m}$ of $\overline{X}_{m}$, the relative error (relerr for short) is defined by

In this subsection, we test the performance of the NNPLS estimator and the RCS estimator for different patterns of fixed basis coefficients. We randomly generated a correlation matrix by the following command:

We took the true matrix $\overline{X}=$ X_bar with dimension n $=500$, rank r $=5$, weight $=5$ and k $=1$. Here, the parameter weight is used to control the relative magnitude difference between the first $k$ largest eigenvalues and the left $r-k$ nonzero eigenvalues. We randomly fixed partial diagonal and off-diagonal entries of $\overline{X}$ and then uniformly sampled the rest entries with i.i.d. Gaussian noise. The noise level, defined by $\|\nu\xi\|_{2}/\|y\|_{2}$ in (2) hereafter, was set to be $10\%$ and the upper bound of the non-fixed diagonal entries was set to be $1$. We further assumed that the rank of the true matrix was known so that for RCS estimator we chose the rank-correction function (23).

In Figure 2, we plot the curves of the relative recovery error and the rank of both the NNPLS estimator (the subfigures on the left) and the RCS estimator (the subfigures on the rigth) for different patterns of fixed entries. Note that both $m$ and $\rho_{m}$ in the rank-correction step (3) depend on the problem of consideration. Thus, we report $m\rho_{m}$ as a whole in the $x$-axis. (Note that for a specific problem, only $\rho_{m}$ is adjustable.) In the captions of subfigures, diag means the number of fixed diagonal entries, and off-diag means the number of fixed off-diagonal entries. For each subfigure on the right side, the initial $\widetilde{X}_{m}$ for the RCS estimator is the point with the smallest recovery error from the corresponding subfigure on the left side.

Figure 2 fully manifests the advantage of the RCS estimator over the NNPLS estimator. It is shown that compared with the NNPLS estimator, the RCS estimator substantially reduces the recovery error and significantly improves the rank consistency. Moreover, the RCS estimator possesses a wide rage of the parameter $\rho_{m}$ to achieve a desired small recovery error and the rank of the true matrix simultaneously. It indicates that whether the resulting solution of a parameter $\rho_{m}$ achieves the true rank can be used to infer the recovery quality. Even if the true rank is unknown in advance, it is still possible to pick out a satisfied solution via monitoring the change of rank in parameter searching. Such advantages are far beyond the reach of the NNPLS estimator.

2 Performance of different rank-correction functions for recovery

3 Performance of different initial NNPLS estimators for recovery

In this subsection, we take the covariance matrix completion for example to test the performance of the RCS estimator with different initial NNPLS estimators $\widetilde{X}_{m}$. We generated the true matrix $\overline{X}$ by the command in Subsection 6.1 with n $=500$, r $=5$, weight $=3$ and k $=1$ except that D = eye(n). The upper bound of the non-fixed diagonal entries was set to be double of the largest absolute value among all the noisy observations of entries together with the fixed entries. We assumed that the rank of the true matrix was known so that we chose the rank-correction function (23).

For each $\rho_{m}$, we first produced the NNPLS estimator, and then use it as the initial point to produce a sequence of RCS estimators with different penalty parameters. Next we choose the RCS estimators that attains the correct rank with the smallest penalty parameter. As can be seen from Figure 2, this choice of the RCS estimator results in the desired small recovery error. The test results are plotted in Figure 4, where the dash curves represent for the NNPLS estimator and the solid curves represent for the chosen RCS estimator. We clearly observe from Figure 4 that, no matter which NNPLS estimator is given to be the initial estimator, the RCS estimator can always substantially improve the recovery quality in terms of both the error and the rank.

4 Performance for different matrix completion problems

In this subsection, we test the performance of the RCS estimator for different matrix completion problems. Figure 2 has revealed that a good choice of the parameter $\rho_{m}$ for the RCS estimator could be the smallest value that attains a stable rank. Therefore, the bisection search method can be used to find such a parameter $\rho_{m}$. This is actually what we benefit from rank consistency. In the following experiments, we apply this strategy to find a suitable $\rho_{m}$ for the RCS estimator.

A natural question then arises: Will multiple rank-correction steps further improve the recovery quality? The answer can be found in Tables 2, 3 and 4 below, which report the experimental results for covariance matrix completion, rectangular matrix completion and density matrix completion, respectively. The reported NNPLS estimator is the one with the smallest recovery error among all different $\rho_{m}$ presuming the true matrix is known. The initial estimator of the first RCS estimator is the NNPLS estimator with a single preset $\rho_{m}=0.4\frac{\eta\|y\|_{2}}{\sqrt{m}}\sqrt{\frac{\log(n_{1}+n_{2})}{mn}}$, where $\eta$ is the noise level. This choice of $\rho_{m}$ follows (9) with $C=0.4$, $\kappa=1$, $\mu=1$ and $\nu$ taken its expected value based on observations. The second (third) RCS estimator takes the first (second) RCS estimator to be the initial estimator. The rank-correction function $F$ is defined by (24), (25) and (26) with $\varepsilon=0.05$ and $\tau=2$.

For the covariance matrix completion problems, we generated the true matrix $\overline{X}$ by the command in Subsection 6.1 with n $=1000$, weight $=2$ and k $=1$ except that D = eye(n). The rank of $\overline{X}$ and the number of fixed diagonal and non-diagonal entries of $\overline{X}$ are reported in the first and the second columns of Table 2, respectively. We sampled partial off-diagonal entries uniformly at random with i.i.d. Gaussian noise at the noise level $10\%$. The upper bound of the non-fixed diagonal entries was set to be double of the largest absolution value among all the noisy observations of entries together with the fixed entries. From Table 2, we see that when the sample ratio is reasonable, a single rank-correction step is fully capable to yield a desired result. However, when the sample ratio is very low, especially if some off-diagonal entries are fixed, one or two further rank-correction steps could still bring some improvement in recovery quality.

For the density matrix completion problems, we generated the true density matrix $\overline{X}$ by the following command:

During the testing, we set n $=1024$, weight $=2$ and k $=1$, and sampled partial Pauli measurements except the trace of $\overline{X}$ uniformly at random with $10\%$ i.i.d. Gaussian noise. Besides this statistical noise, we further added the depolarizing noise, which frequently appears in quantum systems. The strength of the depolarizing noise was set to be $0.01$. This case is labeled as the mixed noise in the last four rows of Table 3. We remark here that the depolarizing noise differs from our assumption on noise since it does not have randomness. One may refer to for details of the quantum depolarizing channel. In , Flammia et al. proposed a two-step method for seeking a feasible solution of low-rank — (1) evaluating an NNPLS estimator by dropping the trace one constraint; (2) normalizing the resulting solution to be of trace one. We tested this method in our experiments, with the NNPLS estimator without trace one constraint chosen to be the one with the smallest recovery error among all that attain the true rank, presuming that the true matrix is known. The two-step results are reported as NNPLS1 and NNPLS2, respectively, in Table 3. Besides the relative recovery error (relerr), we also report the (squared) fidelity, which is a measure of the closeness of two quantum states defined by \big{\|}\widehat{X}_{m}^{1/2}\overline{X}^{1/2}\big{\|}_{*}^{2}. From Table 3, we can see that the RCS estimator is superior to the NNPLS2 estimator in terms of both the fidelity and the relative error.

For the rectangular matrix completion problems, we generated the true matrix $\overline{X}$ by the following command:

We set weight $=2$, k $=1$ and took $\overline{X}=$ X_bar with different dimensions and ranks. Both the uniform sampling scheme and the non-uniform sampling scheme were tested for comparison. For the non-uniform sampling scheme, the probability to sample the first $1/4$ rows and the first $1/4$ columns were $3$ times as much as that of other rows and columns respectively. In other words, the density of sampled entries in the top-left part was $3$ times as much as that in the bottom-left part and the top-right part respectively and $9$ times as much as that in the bottom-right part. We added $10\%$ i.i.d. Gaussian noise to the sampled entries. We also fixed partial entries of $\overline{X}$ uniformly from the rest un-sampled entries. The upper bound of the non-fixed entries was set to be double of the largest absolution value among all the noisy observations of entries together with the fixed entries. What we observe from Table 4 for the rectangular matrix completion is similar to that for the covariance matric completion. Moreover, we can see that the non-uniform sampling scheme greatly weakens the recoverability of the NNPLS estimator in terms of both the recovery error and the rank, especially when the sample ratio is low. Meanwhile, the advantage of the RCS estimators in such cases becomes more remarkable.

Conclusions

In this paper, we proposed a rank-corrected procedure for low-rank matrix completion problems with fixed basis coefficients. This approach can substantially overcome the limitation of the nuclear norm technique for recovering a low-rank matrix. We confirmed the improvement of the rank-correction step in both the reduction of recovery error and the achievement of rank consistency (in the sense of Bach ). Due to the presence of fixed basis coefficients, constraint nondegeneracy plays an important role in our analysis. Extensive numerical experiments show that our approach can significantly improve the recovery performance compared with the nuclear norm penalized least square estimator. As a byproduct, our results also provide a theoretical foundation for the majorized penalty method of Gao and Sun and Gao for structured low-rank matrix optimization problems.

Our proposed rank-correction step also allows additional constraints according to other possible prior information. In order to better fit the under-sampling setting of matrix completion, in the future work, it would be of great interest to extend the asymptotic rank consistency results to the case where the matrix size is allowed to grow. It would also be interesting to extend this approach to deal with other low-rank matrix problems.

Acknowledgements

The authors would like to thank Professor Wotao Yin for his valuable comments on possibly choosing the optimal penalty parameter for recovery error bounds and Dr. Kaifeng Jiang for helpful discussions on efficiently solving the density matrix completion problem.

Appendix Appendix A Spectral operator

where a signed permutation matrix is a real matrix that contains exactly one nonzero entry $1$ or $-1$ in each row and column and elsewhere. From this definition, we see that

Appendix Appendix B Constraint nondegeneracy

Consider the following constrained optimization problem

where $\mathcal{T}_{K}(\widehat{z})$ denotes the tangent cone of $K$ at $\widehat{z}$ and ${\rm lin}(\mathcal{T}_{K}(\widehat{z}))$ denotes the largest linearity space contained in $\mathcal{T}_{K}(\widehat{z})$, i.e., ${\rm lin}(\mathcal{T}_{K}(\widehat{z}))=\mathcal{T}_{K}(\widehat{z})\cap(-\mathcal{T}_{K}(\widehat{z}))$. When the function $\Psi$ is nondifferentiable, we can rewrite the optimization problem (28) equivalently as

From (29) and [67, Theorem 6.41], the constraint nondegeneracy holds at $(\widehat{X},\widehat{t})$ with $\widehat{t}=\Psi(\widehat{X})$ if

By the definition of $\widetilde{\mathcal{A}}$, it is not difficult to verify that this condition is equivalent to

Appendix Appendix C Proofs of Theorems

Let $\Delta_{m}:=\widehat{X}_{m}-\overline{X}$. Using the optimality of $\widehat{X}_{m}$ to the problem (3), we obtain that

where ${\rm row}(X)$ and ${\rm col}(X)$ denote the row space and column space of $X$, respectively. Let $\mathcal{P}_{T}$ and $\mathcal{P}_{T^{\perp}}$ be orthogonal projections onto $T$ and $T^{\perp}$, respectively, given by

Moreover, from the directional derivative of the nuclear norm at $\overline{X}$, (see [75, Theorem 1]), we have

Then, by substituting (33) and (34) into (31), we have

Note that $\text{rank}(\mathcal{P}_{T}(\Delta_{m}))\leq 2r$. Hence, $\|\mathcal{P}_{T}(\Delta_{m})\|_{*}\leq\sqrt{2r}\|\mathcal{P}_{T}(\Delta_{m})\|_{F}\leq\sqrt{2r}\|\Delta_{m}\|_{F}$ and then the desired result (7) follows.

C.2 Proof of Theorem 2

We first show that the sampling operator $\mathcal{R}_{\Omega}$ satisfies some RIP-like property for matrices specified in a certain set with high probability. Similar results can also be found in .

where $\{\epsilon_{1},\ldots,\epsilon_{m}\}$ is an i.i.d. Rademacher sequence, i.e., an i.i.d. sequence of Bernoulli random variables taking the values $1$ and $-1$ with probability $1/2$.

Then, for any given $\gamma>1$, $\tau_{1}\in(0,1)$ and $\tau_{2}\in(0,\tau_{1}/\gamma)$, with probability at least $1-\frac{\exp(-(\tau_{1}-\gamma\tau_{2})^{2}mt^{2}/2)}{1-\exp(-(\gamma^{2}-1)(\tau_{1}-\gamma\tau_{2})^{2}mt^{2}/2)}$,

Proof: The proof is similar to that of [40, Lemma 12]. For any $s,t>0$, $\gamma>1$, $\tau_{1}\in(0,1)$ and $\tau_{2}\in(0,\tau_{1}/\gamma)$, we need to show that the event

occurs with probability less than $\frac{\exp(-(\tau_{1}-\gamma\tau_{2})^{2}mt^{2}/2)}{1-\exp(-(\gamma^{2}-1)(\tau_{1}-\gamma\tau_{2})^{2}mt^{2}/2)}$. We decompose $K(s,t)$ as

For any $a\geq t$, we further define $K(s,t,a):=\{\Delta\in K(s,t)\mid\langle\mathcal{Q}_{\beta}(\Delta),\Delta\rangle\leq a\}.$ Then we get $E\subseteq\bigcup_{k=1}^{\infty}E_{k}$ with

Now we need to estimate the probability of each event $E_{k}$. Define

Since $\|\mathcal{R}_{\beta}(\Delta)\|_{\infty}\leq 1$ for all $\Delta\in K(s,t)$, from Massart’s Hoeffding type concentration inequality [51, Theorem 1.4] for suprema of empirical processes, we have

where the first inequality follows from the symmetrization theorem (e.g., see [73, Lemma 2.3.1] and [6, Theorem 14.3]) and the second inequality follows from the contraction theorem (e.g., see [46, Theorem 4.12] and [6, Theorem 14.4]). Moreover, from (8), we have

Combining (39) and (41) with the definition of $\vartheta_{m}$ in (36), we obtain that

where the second inequality follows from the simple fact $x_{1}x_{2}\leq(x_{1}^{2}+x_{2}^{2})/2$ for any $x_{1},x_{2}\geq 0$. Then, it follows from (38) that

This implies that {\rm Pr}(E_{k})\leq\exp\Big{(}\!\!-\!\frac{1}{2}\,\gamma^{2(k-1)}(\tau_{1}-\gamma\tau_{2})^{2}mt^{2}\Big{)}. Then, since $\gamma>1$, by using $\gamma^{k}\geq 1+k(\gamma-1)$ for any $k\geq 1$, we have

Thus, we complete the proof of Lemma 11. $\Box$

Now we proceed with the proof of Theorem 2. Let $\Delta_{m}:=\widehat{X}_{m}-\overline{X}$. Notice that the equality (35) implies that

This, together with $\|\mathcal{P}_{T}(\Delta_{m})\|_{*}\leq\sqrt{2r}\|\Delta_{m}\|_{F}$, leads to

Let $b_{m}:=\|\mathcal{R}_{\beta}(\Delta_{m})\|_{\infty}\leq 2b$. For any fixed $c>0$, $\gamma>1$, $\tau_{1}\in(0,1)$ and $\tau_{2}\in(0,\tau_{2}/\gamma)$, define $t_{m}:=\sqrt{\frac{2c\log(n_{1}+n_{2})}{(\tau_{1}-\gamma\tau_{2})^{2}m}}$ so that direct calculation yields

Then we separate the discussion into two cases:

Case 1: $\langle\mathcal{Q}_{\beta}(\Delta_{m}),\Delta_{m}\rangle\leq b_{m}^{2}t_{m}$. It follows from (40) that $\|\Delta_{m}\|_{F}^{2}/d_{2}\leq 4b^{2}\mu_{1}t_{m}$.

Case 2: $\langle\mathcal{Q}_{\beta}(\Delta_{m}),\Delta_{m}\rangle>b_{m}^{2}t_{m}$. It follows from (42) that $\Delta_{m}/b_{m}\in K(s_{m},t_{m})$ with s_{m}:=\frac{\kappa}{\kappa-1}\big{(}\sqrt{2}+a_{m}\big{)}\sqrt{r}. Then for any given $\tau_{3}$ satisfying $0<\tau_{3}<1$, we obtain that with probability at least $1-\frac{(n_{1}+n_{2})^{-c}}{1-2^{-(\gamma^{2}-1)c}}$,

where the first inequality follows from (40), the second inequality follows from Lemma 11 and the third inequality follows from Theorem 1. Plugging in $s_{m}$ further leads to

Combing the above two cases together, with $\gamma$, $\tau_{1}$, $\tau_{2}$ and $\tau_{3}$ chosen to be absolute constants, we arrive at an intermediate result that there exist some positive absolute constants $c^{\prime}_{0},c^{\prime}_{1},c^{\prime}_{2}$ and $C^{\prime}_{0}$ such that for any $\kappa>1$, if $\rho_{m}$ is chosen as in Theorem 1, then with probability at least $1-c^{\prime}_{1}(n_{1}+n_{2})^{-c^{\prime}_{2}}$,

Then, there exists a constant $C$ such that for all $t>0$, with probability at least $1\!-\exp(-t)$,

With the help of Lemma 12, we obtain the following result, which is an extension of [44, Lemma 2] and [40, Lemmas 5 & 6] from the standard basis to an arbitrary orthonormal basis. A similar result can also be found in [58, Lemma 6].

Under Assumption 2, there exists a positive constant $C^{\prime}$ (only depending on the $\psi_{1}$ Orlicz norm of $\xi_{k}$) such that for all $t>0$, with probability at least $1-\exp(-t)$,

In particular, when $m\geq\sqrt{d_{2}}\log^{3}(n_{1}+n_{2})/\mu_{2}$, we also have

A good estimation of $\rho_{m}$ can be achieved by choosing $t=c^{\prime}_{2}\log(n_{1}+n_{2})$ in Lemma 13 for an optimal order bound, where $c^{\prime}_{2}$ is the same as that in (43). With this choice, when $m\geq 4(1+c^{\prime}_{2})\sqrt{d_{2}}\log^{2}(d_{2})\log(n_{1}+n_{2})/\mu_{2}$, the first term in the maximum of (44) dominates the second one. Thus, with probability at least $1-(n_{1}+n_{2})^{-c^{\prime}_{2}}$, one can choose

Moreover, since Bernoulli random variables are sub-exponential, Lemma 13 also provides an upper bound of $\vartheta_{m}$ in (45). It is worthwhile to note that after plugging the above estimations of $\rho_{m}$ and $\vartheta_{m}$, the second term in the maximum of (43) is negligible compared with the first term. Therefore, the second term is further dropped for simplicity and thus we complete the proof.

C.3 Proof of Theorem 3

Note that for any $X\in\mathcal{N}_{\delta}(\overline{X})$,

Since $\delta/\sigma_{r}(\overline{X})<1/2$, we further have $\sigma_{r}(X)\geq\sigma_{r}(\overline{X})-\delta>\delta\geq\sigma_{r+1}(X)$. This means

Moreover, $\widehat{F}$ is continuously differentiable over $\mathcal{N}_{\delta}(\overline{X})$. Hence, we can apply the Mean Value Theorem to obtain

where $\widetilde{X}_{t}:=\overline{X}+t(\widetilde{X}-\overline{X})$. Clearly, $\widetilde{X}_{t}\in\mathcal{N}_{\delta}(\overline{X})$ when $t\in$.

Here, “$\circ$” stands for the Hadamard product of matrices. Let $\Delta$ denote the matrix in the bracket of (48). Moreover, let $\Delta_{\chi_{i},\chi_{j}}$ and $\widetilde{H}_{\chi_{i},\chi_{j}}$ denote the submatrices of $\Delta$ and $\widetilde{H}$ with row indices $\chi_{i}$ and column indices $\chi_{j}$, respectively. Then, a direct calculation yields

Note that $\|\widehat{F}^{\prime}(X)(H)\|_{F}=\|\Delta\|_{F}$ and $\|\widetilde{H}\|_{F}=\|H\|_{F}$. By summing up the above inequalities together, we obtain that for any $X\in\mathcal{N}_{\delta}(\overline{X})$,

Now, we proceed with the proof by applying (49) to (47). This leads to

Moreover, using [4, Theorems IV.3.4 & II.3.1], we have

This implies that $\sigma_{r}(\widetilde{X}_{t})-\sigma_{r}(\overline{X})=\delta_{t}\cos\theta$ and $\sigma_{r+1}(\widetilde{X}_{t})=\delta_{t}\sin\theta$ for some $\delta_{t}\leq t\delta$ and $\theta\in[0,2\pi)$. Thus,

Substituting (51) into (50), we obtain that

This, together with (46), completes the proof.

C.4 Proof of Theorem 4

We first prove the following properties of the sample operator $\mathcal{R}_{\Omega}$ and its adjoint $\mathcal{R}_{\Omega}^{*}$.

To prove the convergence in distribution of minimizers, the following theorem of Knight [41, Theorem 1] on epi-convergence in distribution is particularly useful in this regard (see also [32, Proposition 9]).

Let $\{\Phi_{m}\}$ be a sequence of random lower-semicontinuous functions that epi-converges in distribution to $\Phi$. Assume that

$\widehat{x}_{m}$ is an $\varepsilon_{m}$-minimizer of $\Phi_{m}$, i.e., $\Phi_{m}(\widehat{x}_{m})\leq\inf\Phi_{m}(x)+\varepsilon_{m}$, where $\varepsilon_{m}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$;

the function $\Phi$ has a unique minimizer $\overline{x}$.

Then, $\widehat{x}_{m}\stackrel{{\scriptstyle d}}{{\rightarrow}}\overline{x}$. In addition, if $\Phi$ is a deterministic function, then $\widehat{x}_{m}\stackrel{{\scriptstyle p}}{{\rightarrow}}\overline{x}$.

It is know from that $\widehat{x}_{m}$ is guaranteed to be $O_{p}(1)$ when all $\Phi_{m}$ are convex functions and $\Phi$ has a unique minimizer. For more details on epi-convergence in distribution, one may refer to King and Wets , Geyer , Pflug and Knight . As Lemma 15 is only applicable to unconstrained optimization problems, constrained optimization problems need to be equivalently converted to unconstrained ones using the indicator function of feasible set. This leads to the issue of epi-convergence in distribution of the sum of two sequences of random functions; see, e.g., Pflug [60, Lemma 1].

Now we proceed with the proof of Theorem 4. Let $\Phi_{m}$ denote the objective function of (3) and $\mathcal{F}$ denote the feasible set. Then, the problem (3) can be concisely written as

C.5 Proof of Theorem 5

Theorem 4 actually implies that $\widehat{X}_{m}$ has a higher rank than $\overline{X}$ with probability converging to $1$ if $\rho_{m}\rightarrow 0$, due to the straightforward result:

If $X_{m}\stackrel{{\scriptstyle p}}{{\rightarrow}}\overline{X}$, then \lim\limits_{m\rightarrow\infty}{\rm Pr}\big{(}{\rm rank}(X_{m})\geq{\rm rank}(\overline{X})\big{)}=1.

Proof: It follows from the Lipschitz continuity of singular values that

Now we take a look at the local property for the rank function for the perturbation.

This implies that ${\rm rank}(\overline{X}+\rho\Delta)>{\rm rank}(\overline{X})$. $\Box$

Define $\widehat{\Delta}_{m}:=\rho_{m}^{-1}(\widehat{X}_{m}-\overline{X})$. To guarantee the efficiency of the nuclear semi-norm on encouraging a low-rank solution, the parameter $\rho_{m}$ should not decay too fast. Then, for a slow decay on $\rho_{m}$, we can establish the following result.

If $\rho_{m}\rightarrow 0$ and $\sqrt{m}\rho_{m}\rightarrow\infty$, then $\widehat{\Delta}_{m}\stackrel{{\scriptstyle p}}{{\rightarrow}}\widehat{\Delta}$, where $\widehat{\Delta}$ is the unique optimal solution to the following convex optimization problem

Proof: Take a variable transformation $\Delta:=\rho_{m}^{-1}(X-\overline{X})$ in the optimization problem (3). Then one can easily see that $\widehat{\Delta}_{m}$ is the optimal solution to

Then, under Assumptions 3 and 4, according to Lemma 14, we obtain that $\Phi_{m}$ converges pointwise in probability to $\Phi$. Together with the convexity of $\mathcal{K}$, we know that $\mathcal{F}_{m}$ converges in the sense of Painlevé-Kuratowski to the tangent cone $\mathcal{T}_{\mathcal{K}}(\overline{X})$ (see ), taking the form

Since epi-convergence of functions corresponds to set convergence of their epigraphs , we obtain that $\delta_{\mathcal{F}_{m}}$ epi-converges to $\delta_{\mathcal{T}_{\mathcal{K}}(\overline{X})}$. Then, by using the same argument as in the proof of Theorem 4, we obtain that $\Phi_{m}+\delta_{\mathcal{F}_{m}}$ epi-converges in distribution to $\Phi+\delta_{\mathcal{T}_{\mathcal{K}}(\overline{X})}$. In addition, the optimal solution to (52) is unique due to the strong convexity of $\Phi$ over the feasible set $\mathcal{K}$. Then, applying Lemma 15 on the epi-convergence in distribution leads to the desired result. $\Box$

Note that $\mathcal{R}_{\beta^{+}}(\widehat{\Delta})\leq 0$ and $\mathcal{R}_{\beta^{-}}(\widehat{\Delta})\geq 0$ implies that $\mathcal{Q}_{\beta}^{\dagger}\mathcal{Q}_{\beta}(\widehat{\Delta})=\mathcal{P}_{\beta}(\widehat{\Delta})$. Moreover, $\mathcal{Q}_{\beta}^{\dagger}\mathcal{R}_{\alpha}^{*}(\widehat{\eta}^{0})=\mathcal{Q}_{\beta}^{\dagger}\mathcal{R}_{\beta^{+}}^{*}(\widehat{\eta}^{1})=\mathcal{Q}_{\beta}^{\dagger}\mathcal{R}_{\beta^{-}}^{*}(\widehat{\eta}^{2})=0$. Then, we apply the operator $\mathcal{Q}_{\beta}^{\dagger}$ to the first equation of (55) and then obtain

Then, by further applying the operator $\mathcal{Q}_{\beta}^{\dagger}$ to the above equation, together with (56) in Lemma 19 and (57), we obtain that

C.6 Proof of Theorem 6

The proof of Theorem 6 is similar to the proof of Theorem 5. Define $\widehat{\Delta}_{m}:=\rho_{m}^{-1}(\widehat{X}_{m}-\overline{X})$.

If $\rho_{m}\rightarrow 0$ and $\sqrt{m}\rho_{m}\rightarrow\infty$, then $\widehat{\Delta}_{m}\stackrel{{\scriptstyle p}}{{\rightarrow}}\widehat{\Delta}$, where $\widehat{\Delta}$ is the unique optimal solution to the following convex optimization problem

Proof: It is easy to verify that $\widehat{\Delta}_{m}$ is the optimal solution to

C.7 Proof of Theorem 7

C.8 Proof of Theorem 9

We first prove for the constraint nondegeneracy.

For the matrix completion problems of Classes I and II, the constraint nondegeneracy (16) holds at $\overline{X}$.

Proof: For the real covariance matrix case, the proof is given in [61, Lemma 3.3] and [62, Proposition 2.1]. For the complex covariance matrix case, one can use the similar arguments to prove the result.

This means that the constraint nondegeneracy (16) holds. $\Box$

Since $\text{rank}(\overline{X})=r$, with the choice (22) of $F$, we have that for any $i\in\pi$,

By taking the trace on both sides, we obtain that $\widehat{\Lambda}=\frac{1}{r}\text{Tr}(F(\overline{X}))I_{n-r}.$ Since $\overline{X}$ is a density matrix of rank $r$, with the choice (22) of $F$, we have that

References