Inference and Uncertainty Quantification for Noisy Matrix Completion

Introduction

Low-rank matrix completion is concerned with recovering a low-rank matrix, when only a small fraction of its entries are revealed to us [Sre04, CR09, KMO10a]. Tackling this problem in large-scale applications is computationally challenging, due to the intrinsic nonconvexity incurred by the low-rank structure. To further complicate matters, another inevitable challenge stems from the imperfectness of data acquisition mechanisms, wherein the acquired samples are contaminated by a certain amount of noise.

Fortunately, if the entries of the unknown matrix are sufficiently de-localized and randomly revealed, this problem may not be as hard as it seems. Substantial progress has been made over the past several years in designing computationally tractable algorithms — including both convex and nonconvex approaches — that allow to fill in unseen entries faithfully given only partial noisy samples [CP10, NW12, KLT11, KMO10b, CW15, MWCC17, CCF+19]. Nevertheless, modern decision making would often require one step further. It not merely anticipates a faithful estimate, but also seeks to quantify the uncertainty or “confidence” of the provided estimate, ideally in a reasonably accurate fashion. For instance, given an estimate returned by the convex approach, how to use it to compute a short interval that is likely to contain a missing entry?

Conducting effective uncertainty quantification for noisy matrix completion is, however, far from straightforward. For the most part, the state-of-the-art matrix completion algorithms require solving highly complex optimization problems, which often do not admit closed-form solutions. Of necessity, it is generally very challenging to pin down the distributions of the estimates returned by these algorithms. The lack of distributional characterizations presents a major roadblock to performing valid, yet efficient, statistical inference on the unknown matrix of interest.

It is worth noting that a number of recent papers have been dedicated to inference and uncertainty quantification for various high-dimensional problems in high-dimensional statistics, including Lasso [ZZ14, vdGBRD14, JM14a], generalized linear models [vdGBRD14, NL17, BFL+18], graphical models [JVDG15, RSZZ15, MLL17]), amongst others. Very little work, however, has looked into noisy matrix completion along this direction. While non-asymptotic statistical guarantees for noisy matrix completion have been derived in prior theory, most, if not all, of the estimation error bounds are supplied only at an order-wise level. Such order-wise error bounds either lose a significant factor relative to the optimal guarantees, or come with an unspecified (but often enormous) pre-constant. Viewed in this light, a confidence region constructed directly based on such results is bound to be overly conservative, resulting in substantial over-coverage.

2 A glimpse of our contributions

This paper takes a substantial step towards efficient inference and uncertainty quantification for noisy matrix completion. Specifically, we develop a simple procedure to compensate for the bias of the commonly used convex and nonconvex estimators. The resulting de-biased estimators admit nearly accurate non-asymptotic distributional guarantees. Such distributional characterizations in turn allow us to reason about the uncertainty of the obtained estimates vis-à-vis the unknown matrix. While details of our main findings are postponed to Section 3, we would like to immediately single out a few important merits of the proposed inferential procedures and theory:

Our results enable two types of uncertainty assessment, namely, we can construct (i) confidence intervals for each entry — either observed or missing — of the unknown matrix; (ii) confidence regions for the low-rank factors of interest (modulo some unavoidable global ambiguity).

Despite the complicated statistical dependency, our procedure and theory do not rely on sample splitting, thus avoiding the unnecessary widening of confidence intervals / regions due to insufficient data usage.

The confidence intervals / regions constructed based on the proposed procedures are, in some sense, optimal.

We present a unified approach that accommodates both convex and nonconvex estimators seamlessly.

As a byproduct, we characterize the Euclidean estimation errors of the proposed de-biased estimators. Such error bounds are sharp and match an oracle lower bound precisely (including the pre-constant). To the best of our knowledge, this is the first theory that demonstrates that a computationally feasible algorithm can achieve the statistical limit including the pre-constant.

All of this is built upon the intimate link between convex and nonconvex estimators [CCF+19], as well as the recent advances in analyzing the stability of nonconvex optimization against random noise [MWCC17].

Models and notation

To cast the noisy matrix completion problem in concrete statistical settings, we adopt a model commonly studied in the literature [CR09]. We also introduce some useful notation.

Observation models.

What we observe is a random subset of noisy entries of $\bm{M}^{\star}$; more specifically, we observe

Incoherence conditions.

Clearly, not all matrices can be reliably estimated from a highly incomplete set of measurements. To address this issue, we impose a standard incoherence condition [CR09, Che15] on the singular subspaces of $\bm{M}^{\star}$ (i.e. $\bm{U}^{\star}$ and $\bm{V}^{\star}$):

Asymptotic notation.

Here, $f(n)\lesssim h(n)$ (or $f(n)=O(h(n))$) means $|f(n)|\leq c_{1}|h(n)|$ for some constant $c_{1}>0$, $f(n)\gtrsim h(n)$ means $|f(n)|\geq c_{2}|h(n)|$ for some constant $c_{2}>0$, $f(n)\asymp h(n)$ means $c_{2}|h(n)|\leq|f(n)|\leq c_{1}|h(n)|$ for some constants $c_{1},c_{2}>0$, and $f(n)=o(h(n))$ means $\lim_{n\rightarrow\infty}f(n)/h(n)=0$. We write $f(n)\ll h(n)$ to indicate that $|f(n)|\leq c_{1}|h(n)|$ for some small constant $c_{1}>0$ (much smaller than 1), and use $f(n)\gg h(n)$ to indicate that $|f(n)|\geq c_{2}|h(n)|$ for some large constant $c_{2}>0$ (much larger than 1).

Inferential procedures and main results

The proposed inferential procedure lays its basis on two of the most popular estimation paradigms — convex relaxation and nonconvex optimization — designed for noisy matrix completion. Recognizing the complicated bias of these two highly nonlinear estimators, we shall first illustrate how to perform bias correction, followed by a distributional theory that establishes the near-Gaussianity and optimality of the proposed de-biased estimators.

We first review in passing two tractable estimation algorithms that are arguably the most widely used in practice. They serve as the starting point for us to design inferential procedures for noisy low-rank matrix completion. The readers familiar with this literature can proceed directly to Section 3.2.

Recall that the rank function $\mathsf{rank}(\cdot)$ is highly nonconvex, which often prevents us from computing a rank-constrained estimator in polynomial time. For the sake of computational feasibility, prior works suggest relaxing the rank function into its convex surrogate [Faz02, RFP10]; for example, one can consider the following penalized least-squares convex program

or using our notation $\mathcal{P}_{\Omega}$,

Here, $\|\cdot\|_{*}$ is the nuclear norm (the sum of singular values, which is a convex surrogate of the rank function), and $\lambda>0$ is some regularization parameter. Under mild conditions, the solution to the convex program (3.1) provably attains near-optimal estimation accuracy (in an order-wise sense), provided that a proper regularization parameter $\lambda$ is adopted [CCF+19].

Nonconvex optimization.

Intimate connections between convex and nonconvex estimates.

Denote by $\bm{Z}^{\mathsf{cvx}}$ any minimizer of the convex program (3.1), and denote by $(\bm{X}^{\mathsf{ncvx}},\bm{Y}^{\mathsf{ncvx}})$ the estimate returned by Algorithm 1 aimed at solving (3.4). As was recently shown in [CCF+19], when the regularization parameter $\lambda$ is properly chosen, these two estimates obey (see (A.12) in Appendix A.2 for a precise statement)

2 Constructing de-biased estimators

We are now well equipped to describe how to construct new estimators based on the convex estimate $\bm{Z}^{\mathsf{cvx}}$ and the nonconvex estimate $(\bm{X}^{\mathsf{ncvx}},\bm{Y}^{\mathsf{ncvx}})$, so as to enable statistical inference. Motivated by the proximity of the convex and nonconvex estimates and for the sake of conciseness, we shall abuse notation by using the shorthand $\bm{Z},\bm{X},\bm{Y}$ for both convex and nonconvex estimates; see Table 1 and Appendix B for precise definitions. This allows us to unify the presentation for both convex and nonconvex estimators.

Given that both (3.1) and (3.4) are regularized least-squares problems, they behave effectively like shrinkage estimators, indicating that the provided estimates necessarily suffer from non-negligible bias. In order to enable desired statistical inference, it is natural to first correct the estimation bias.

A natural de-biasing strategy that immediately comes to mind is the following simple linear transformation (recall the notation in Table 1):

To remedy this issue, we propose to further project $\bm{Z}^{0}$ onto the set of rank-$r$ matrices, leading to the following de-biased estimator

The estimator (3.7) can be viewed as performing one iteration of singular value projection (SVP) [MJD09, DC18] on the current estimate $\bm{Z}$.

The estimator (3.7) also bears a similarity to the de-biased estimator proposed by [Xia18] for low-rank trace regression; the disparity between them shall be discussed in Section 4.

An equivalent form: a de-shrunken estimator for the low-rank factors.

It turns out that the de-biased estimator (3.7) admits another almost equivalent representation that offers further insights. Specifically, we consider the following de-shrunken estimator for the low-rank factors

where we recall the definition of $\bm{X}$ and $\bm{Y}$ in Table 1. To develop some intuition regarding why this is called a de-shrunken estimator, let us look at a simple scenario where $\bm{U}\bm{\Sigma}\bm{V}^{\top}$ is the SVD of $\bm{X}\bm{Y}^{\top}$ and $\bm{X}=\bm{U}\bm{\Sigma}^{1/2}$, $\bm{Y}=\bm{V}\bm{\Sigma}^{1/2}$. It is then self-evident that

In words, $\bm{X}^{\mathsf{d}}$ and $\bm{Y}^{\mathsf{d}}$ are obtained by de-shrinking the spectrum of $\bm{X}$ and $\bm{Y}$ properly.

As we shall formalize in Section 5.1, the de-shrunken estimator (3.8) for the low-rank factors is nearly equivalent to the de-biased estimator (3.7) for the whole matrix, in the sense that

3 Main results: distributional guarantees

The proposed estimators admit tractable distributional characterizations in the large-$n$ regime, which facilitates the construction of confidence regions for many quantities of interest. In particular, this paper centers around two types of inferential problems:

Each entry of the matrix $\bm{M}^{\star}$: the entry can be either missing (i.e. predicting an unseen entry) or observed (i.e. de-noising an observed entry). For example, in the Netflix challenge, one would like to infer a user’s preference about any movie, given partially revealed ratings [CR09]. Mathematically, this seeks to determine the distribution of

Suppose that the sample size and the noise obey

A more complete version can be found in Theorem 5.

Another interesting feature — which we shall make precise in the proof of this theorem — is that: for any given $1\leq i,j\leq n$, the two random vectors $\bm{Z}_{\bm{X}}^{\top}\bm{e}_{i}$ and $\bm{Z}_{\bm{Y}}^{\top}\bm{e}_{j}$ are nearly statistically independent. This is crucial for deriving inferential guarantees for the entries of the matrix.

Theorem 1 is a non-asymptotic result. In words, Theorem 1 decomposes the estimation error $\bm{X}^{\mathsf{d}}\bm{H}^{\mathsf{d}}-\bm{X}^{\star}$ (resp. $\bm{Y}^{\mathsf{d}}\bm{H}^{\mathsf{d}}-\bm{Y}^{\star}$) into a Gaussian component $\bm{Z}_{\bm{X}}$ (resp. $\bm{Z}_{\bm{Y}}$) and a residual term $\bm{\Psi}_{\bm{X}}$ (resp. $\bm{\Psi}_{\bm{Y}}$). If the sample size is sufficiently large and the noise size is sufficiently small, then the residual terms are much smaller in size compared to $\bm{Z}_{\bm{X}}$ and $\bm{Z}_{\bm{Y}}$. To see this, it is helpful to leverage the Gaussianity (3.15a) and compute that: for each $1\leq j\leq n$, the $j$th row of $\bm{Z}_{\bm{X}}$ obeys

in other words, the typical size of the $j$th row of $\bm{Z}_{\bm{X}}$ is no smaller than the order of $\sigma\sqrt{r/(p\sigma_{\max})}$. In comparison, the size of each row of $\bm{\Psi}_{\bm{X}}$ (see (3.16)) is much smaller than $\sigma\sqrt{r/(p\sigma_{\max})}$ (and hence smaller than the size of the corresponding row of $\bm{Z}_{\bm{X}}$) with high probability, provided that (3.13) is satisfied.

Equipped with the above master decompositions of the low-rank factors and Remark 4, we are ready to present a similar decomposition for the entry $M_{ij}^{\mathsf{d}}-M_{ij}^{\star}$.

For each $1\leq i,j\leq n$, define the variance $v_{ij}^{\star}$ as

where $\bm{U}_{i,\cdot}^{\star}$ (resp. $\bm{V}_{j,\cdot}^{\star}$) denotes the $i$th (resp. $j$th) row of $\bm{U}^{\star}$ (resp. $\bm{V}^{\star}$). Suppose that

where $g_{ij}\sim\mathcal{N}(0,v_{ij}^{\star})$ and the residual obeys $|\Delta_{ij}|=o(\sqrt{v_{ij}^{\star}})$ with probability exceeding $1-O(n^{-3})$.

In the symmetric case where the noise $\bm{E}$, the truth $\bm{M}^{\star}$, and the sampling pattern are all symmetric (i.e. \mathcal{P}_{\Omega}(\bm{E})=\big{(}\mathcal{P}_{\Omega}(\bm{E})\big{)}^{\top} and $\bm{M}^{\star}=\bm{M}^{\star\top}$), the variance $v_{ii}^{\star}$ (cf. (3.17)) for the diagonal entries has a different formula; more specifically, it is straightforward to extend our theory to show that

This additional multiplicative factor of 2 arises since $\bm{Z}_{\bm{X}}^{\top}\bm{e}_{i}$ and $\bm{Z}_{\bm{Y}}^{\top}\bm{e}_{i}$ are identical (and hence not independent) in this symmetric case. The variance formula for any $v_{ij}^{\star}$ ($i\neq j$) remains unchanged.

Several remarks are in order. To begin with, we develop some intuition regarding where the variance $v_{ij}^{\star}$ comes from. By virtue of Theorem 1, one has the following Gaussian approximation

Assuming that the first-order expansion is reasonably tight, one has

According to Remark 4, $\bm{Z}_{\bm{X}}^{\top}\bm{e}_{i}$ and $\bm{Z}_{\bm{Y}}^{\top}\bm{e}_{j}$ are nearly independent. It is thus straightforward to compute the variance of (3.20) as

Here, (i) relies on (3.20) and the near independence between $\bm{Z}_{\bm{X}}^{\top}\bm{e}_{i}$ and $\bm{Z}_{\bm{Y}}^{\top}\bm{e}_{j}$; (ii) uses the variance formula in Theorem 1; (iii) arises from the definitions of $\bm{X}^{\star}$ and $\bm{Y}^{\star}$ (cf. (2.2)). This computation explains (heuristically) the variance formula $v_{ij}^{\star}$.

Given that Theorem 2 reveals the tightness of Gaussian approximation under conditions (3.18), it in turn allows us to construct nearly accurate confidence intervals for each matrix entry $M_{ij}^{\star}$. This is formally summarized in the following corollary, the proof of which is deferred to Appendix F. Here and throughout, we use $[a\pm b]$ to denote the interval $[a-b,a+b]$.

Denote by $\Phi(t)$ the CDF of a standard Gaussian random variable and by $\Phi^{-1}(\cdot)$ its inverse function. Let

be the empirical estimate of the theoretical variance $v_{ij}^{\star}$. Then one has

In words, Corollary 1 tells us that for any fixed significance level $0<\alpha<1$, the interval

is a nearly accurate two-sided $(1-\alpha)$ confidence interval of $M_{ij}^{\star}$.

In addition, we remark that when $\|\bm{U}_{i,\cdot}^{\star}\|_{2}=\|\bm{V}_{j,\cdot}^{\star}\|_{2}=0$ (and hence $V_{ij}^{\star}=0$), the above Gaussian approximation is completely off. In this case, one can still leverage Theorem 1 to show that

Last but not least, the careful readers might wonder how to interpret our conditions on the sample complexity and the signal-to-noise ratio. Take the case with $r,\mu,\kappa=O(1)$ for example: our conditions read

4 Lower bounds and optimality for inference

It is natural to ask how well our inferential procedures perform compared to other algorithms. Encouragingly, the de-biased estimator is optimal in some sense; for instance, it nearly attains the minimum covariance among all unbiased estimators. To formalize this claim, we shall

Quantify the performance of two ideal estimators with the assistance of an oracle;

Demonstrate that the performance of our de-biased estimators is arbitrarily close to that of the ideal estimators.

In what follows, we denote by $\bm{X}_{i,\cdot}^{\star}$ (resp. $\bm{Y}_{i,\cdot}^{\star}$) the $i$th row of $\bm{X}^{\star}$ (resp. $\bm{Y}^{\star}$).

Suppose that there is an oracle informing us of $\bm{Y}^{\star}$, and that we observe the same set of data as in (2.3). Under such an idealistic setting and for any given $1\leq i\leq n$, the following least-squares estimator achieves the minimum covariance among all unbiased estimators for the $i$th row $\bm{X}_{i,\cdot}^{\star}$ of $\bm{X}^{\star}$ (see e.g. [Sha03, Theorem 3.7])

In other words, for any unbiased estimator $\bm{u}$ of $\bm{X}_{i,\cdot}^{\star}$ (conditional on $\Omega$), one has

where \mathsf{Cov}\big{(}\bm{X}_{i,\cdot}^{\mathsf{ideal}}\,\big{|}\,\Omega\big{)} is precisely the Cramér-Rao lower bound (conditional on $\Omega$) under this ideal setting. As it turns out, with high probability, this lower bound concentrates around $\sigma^{2}(\bm{\Sigma}^{\star})^{-1}/p$, as stated in the following lemma. The proof is postponed to Appendix H.1.

Fix an arbitrarily small constant $\varepsilon>0$. Suppose that $n^{2}p\geq C_{0}\varepsilon^{-2}\kappa^{4}\mu rn$ for some sufficiently large constant $C_{0}>0$ independent of $n$. Then with probability at least $1-O(n^{-10})$, one has

Given that $\varepsilon$ can be an arbitrarily small constant, Lemma 1 uncovers that the covariance of the de-shrunken estimator $\bm{X}_{i,\cdot}^{\mathsf{d}}$ (cf. Theorem 1) matches that of the ideal estimator $\bm{X}_{i,\cdot}^{\mathsf{ideal}}$, thus achieving the Cramér-Rao lower bound with high probability. The same conclusion applies to $\bm{Y}_{j,\cdot}^{\mathsf{d}}$ as well.

Suppose that there is another oracle informing us of $\{\bm{X}_{k,\cdot}^{\star}\}_{k:k\neq i}$ and $\{\bm{Y}_{k,\cdot}^{\star}\}_{k:k\neq j}$; that is, everything about $\bm{X}^{\star}$ except $\bm{X}_{i,\cdot}^{\star}$ and everything about $\bm{Y}^{\star}$ except $\bm{Y}_{j,\cdot}^{\star}$. In addition, we observe the same set of data as in (2.3), except that we do not get to see $M_{ij}$.The exclusion of $M_{ij}$ is merely for ease of presentation. One can consider the model where all $M_{ij}$ with $(i,j)\in\Omega$ are observed with a slightly more complicated argument. Under this idealistic model, the Cramér-Rao lower bound [Sha03, Theorem 3.3] for estimating $M_{ij}^{\star}=\bm{X}_{i,\cdot}^{\star}(\bm{Y}_{j,\cdot}^{\star})^{\top}$ can be computed as

This means that any unbiased estimator of $M_{ij}^{\star}$ must have variance no smaller than $\mathsf{CRLB}(M_{ij}^{\star}\mid\Omega)$. This quantity admits a much simpler lower bound as follows, whose proof can be found in Appendix H.2.

Fix an arbitrarily small constant $\varepsilon>0$. Suppose that $n^{2}p\geq C_{0}\varepsilon^{-2}\kappa^{4}\mu rn\log n$ for some sufficiently large constant $C_{0}>0$ independent of $n$. Then with probability at least $1-O(n^{-10})$,

where $v_{ij}^{\star}$ is defined in Theorem 2.

Similar to Lemma 1, Lemma 2 reveals that the variance of our de-biased estimator $M_{ij}^{\mathsf{d}}$ (cf. Theorem 2) — which certainly does not have access to the side information provided by the oracle — is arbitrarily close to the Cramér-Rao lower bound aided by an oracle.

All in all, the above lower bounds demonstrate that the degrees of uncertainty underlying our de-shrunken low-rank factors and de-biased matrix are, in some sense, statistically minimal.

5 Back to estimation: the de-biased estimator is optimal

While the emphasis of the current paper is on inference, we would nevertheless like to single out an important consequence that informs the estimation step. To be specific, the decompositions and distributional guarantees derived in Theorem 1 and Theorem 2 allow us to track the estimation accuracy of $\bm{M}^{\mathsf{d}}$, as stated in the following theorem. The proof of this result is postponed to Appendix G.

In stark contrast to prior statistical estimation guarantees (e.g. [CP10, NW12, KLT11, CCF+19]), Theorem 3 pins down the estimation error of the proposed de-biased estimator in a sharp manner (namely, even the pre-constant is fully determined). Encouragingly, there is a sense in which the proposed de-biased estimator achieves the best possible statistical estimation accuracy, as revealed by the following result.

Fix an arbitrarily small constant $\varepsilon>0$. Suppose that $n^{2}p\gtrsim\mu rn\log^{2}n$, and that $r=o(n)$. Then with probability exceeding $1-O(n^{-10})$, any unbiased estimator $\widehat{\bm{M}}$ of $\bm{M}^{\star}$ obeys

Intuitively, the term $2nr$ reflects approximately the underlying degrees of freedom in the true subspace $T^{\star}$ of interest (i.e. the tangent space of the rank-$r$ matrices at the truth $\bm{M}^{\star}$), whereas the factor $1/p$ captures the effect due to sub-sampling. This result has already been established in [CP10, Section III.B] (together with [CR09, Theorem 4.1]). We thus omit the proof for conciseness. The key idea is to consider an oracle informing us of the true tangent space $T^{\star}$. ∎

The implication of the above two theorems is remarkable: the de-biasing step not merely facilitates uncertainty assessment, but also proves crucial in minimizing the estimation errors. It achieves optimal statistical efficiency in terms of both the rate and the pre-constant. As far as we know, this is the first theory about a polynomial time algorithm that matches the statistical limit in terms of the pre-constant. This intriguing finding is further corroborated by numerical experiments; see Section 3.6 for details (in particular, Figure 3).

6 Numerical experiments

We conduct numerical experiments on synthetic data to verify the distributional characterizations provided in Theorem 1 and Theorem 2. Note that our main results hold for the de-biased estimators built upon $\bm{Z}^{\mathsf{cvx}}$ and $\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}$. As we will formalize shortly in Section 5.1, these two de-biased estimators are extremely close to each other; see also Figure 4 for experimental evidence. Therefore, in order to save space, we use the de-biased estimator built upon the convex estimate $\bm{Z}^{\mathsf{cvx}}$ throughout the experiments.

We begin by checking the validity of Theorem 1. Suppose that one is interested in estimating the inner product $\bm{e}_{i}^{\top}\bm{X}^{\star}\bm{X}^{\star\top}\bm{e}_{j}$ between $\bm{X}^{\star\top}\bm{e}_{i}$ and $\bm{X}^{\star\top}\bm{e}_{j}$ ($i\neq j$). In the Netflix challenge, this might correspond to the similarity between the $i$th user and the $j$th one. As a straightforward consequence of Theorem 1, the normalized estimation error

is extremely close to a standard Gaussian random variable. Here, similar to (3.22), we let

be the empirical estimate of the theoretically predicted variance $\sigma^{2}(\|\bm{U}_{i,\cdot}^{\star}\|_{2}^{2}+\|\bm{U}_{j,\cdot}^{\star}\|_{2}^{2})/p$. As a result, a $95\%$ confidence interval of $\bm{e}_{i}^{\top}\bm{X}^{\star}\bm{X}^{\star\top}\bm{e}_{j}$ would be $[\bm{e}_{i}^{\top}\bm{X}^{\mathsf{d}}\bm{X}^{\mathsf{d}\top}\bm{e}_{j}\pm 1.96\sqrt{\rho_{ij}}]$. For each $(i,j)$, we define $\widehat{\mathsf{Cov}}_{\mathsf{L},(i,j)}$ to be the empirical coverage rate of $\bm{e}_{i}^{\top}\bm{X}^{\star}\bm{X}^{\star\top}\bm{e}_{j}$ over $200$ Monte Carlo simulations. Correspondingly, denote by $\mathsf{Mean}(\widehat{\mathsf{Cov}}_{\mathsf{L}})$ (resp. $\mathsf{Std}(\widehat{\mathsf{Cov}}_{\mathsf{L}})$) the average (resp. the standard deviation) of $\widehat{\mathsf{Cov}}_{\mathsf{L},(i,j)}$ over indices $1\leq i<j\leq n$. Table 2 collects the simulation results for different values of $(r,p,\sigma)$. As can be seen, the reported empirical coverage rates are reasonably close to the nominal level $95\%$. In addition, Figure 1 depicts the Q-Q (quantile-quantile) plots of $T_{12},T_{13}$ and $T_{14}$ vs. the standard Gaussian random variable over 200 Monte Carlo simulations for $r=5$, $p=0.4$ and $\sigma=10^{-3}$. It is clearly seen that all of these are well approximated by a standard Gaussian random variable.

Next, we turn to Theorem 2, namely the distributional guarantee for the entries of the matrix. Denote

where $v_{ij}$ is the empirical variance defined in (3.22). In view of the $95\%$ confidence interval predicted by Corollary 1, and similar to what have done for the low-rank components, for each $(i,j)$, we define $\widehat{\mathsf{Cov}}_{\mathsf{E},(i,j)}$ to be the empirical coverage rate of $M^{\star}_{ij}$ over $200$ Monte Carlo simulations. Correspondingly, denote by $\mathsf{Mean}(\widehat{\mathsf{Cov}}_{\mathsf{E}})$ (resp. $\mathsf{Std}(\widehat{\mathsf{Cov}}_{\mathsf{E}})$) the average (resp. the standard deviation) of $\widehat{\mathsf{Cov}}_{\mathsf{E},(i,j)}$ over indices $1\leq i,j\leq n$. As before, Table 3 gathers the empirical coverage rates for $M^{\star}_{ij}$ and Figure 2 displays the Q-Q (quantile-quantile) plots of $S_{11},S_{12}$ and $S_{13}$ vs. the standard Gaussian random variable over 200 Monte Carlo trials for $r=5$, $p=0.4$ and $\sigma=10^{-3}$. It is evident that the distribution of $S_{ij}$ matches that of $\mathcal{N}(0,1)$ reasonably well.

We conclude this section with experiments on real data. Similar to [CP10], we use the daily temperature data [NCD19] for 1400 stations across the world in 2018, which results in a $1400\times 365$ data matrix. Inspection on the singular values reveals that the data matrix is nearly low-rank. We vary the observation probability $p$ from $0.5$ to $0.9$ and randomly subsample the data accordingly. Based on the observed temperatures, we then apply the proposed methodology to obtain $95\%$ confidence intervals for all the entries. Table 4 reports the empirical coverage probabilities, the average length of the confidence intervals as well as the estimation error of both $\bm{Z}^{\mathsf{cvx}}$ and $\bm{M}^{\mathsf{d}}$ over 20 independent experiments. It can be seen that the average coverage probabilities are reasonably close to $95\%$ and the confidence intervals are also quite short. In addition, the estimation error of $\bm{M}^{\mathsf{d}}$ is smaller than that of $\bm{Z}^{\mathsf{cvx}}$, which corroborates our theoretical prediction. The discrepancy between the nominal coverage probability and the actual one might arise from the facts that (1) the underlying true temperature matrix is only approximately low-rank, and (2) the noise in the temperature might not be independent.

7 A bit of intuition

We pause to develop some intuition behind the distributional guarantees for the proposed estimators. Bearing in mind the intimate link between convex and nonconvex optimization (cf. (3.5)), it suffices to concentrate on the nonconvex problem (3.4). For the sake of clarity, we further restrict attention to the rank-1 positive semidefinite case where $\bm{M}^{\star}=\bm{x}^{\star}\bm{x}^{\star\top}$ and set $\lambda=0$, where one can focus on

Any optimizer $\widehat{\bm{x}}$ of (3.34) would necessarily satisfy the first-order optimality condition

We shall also assume that $\widehat{\bm{x}}$ is a reasonably reliable estimate obeying $\widehat{\bm{x}}\approx\bm{x}^{\star}$.

We begin with the no-missing-data case (i.e. $p=1$), which already conveys the key insight. The condition (3.35) simplifies to

which, through a little manipulation, leads to an equivalent decomposition:

Then: (1) the third term of (3.37), which can be viewed as a second-order term (in the sense that it is a quadratic term of $\widehat{\bm{x}}-\bm{x}^{\star}$), becomes vanishingly small when $\widehat{\bm{x}}\approx\bm{x}^{\star}$; (2) while the second term of (3.37) looks like a first-order term, it is natural to conjecture that $\widehat{\bm{x}}-\bm{x}^{\star}$ is sufficiently random and hence $(\widehat{\bm{x}}-\bm{x}^{\star})^{\top}\bm{x}^{\star}\ll\|\widehat{\bm{x}}-\bm{x}^{\star}\|_{2}\|\bm{x}^{\star}\|_{2}$ (i.e. the estimation error is not aligned with $\bm{x}^{\star}$), meaning that this term is also expected to be negligible compared to a typical first-order term (e.g. the term on the left-hand side of 3.37). In summary, these non-rigorous arguments suggest that

If one can be convinced that $\bm{E}$ and $\widehat{\bm{x}}$ are only weakly dependent, then this means

Returning to the missing data scenario with $p<1$, everything is based on the following approximation

this is certainly expected — using standard concentration arguments — if we “pretend” that $\mathcal{P}_{\Omega}$ and $\widehat{\bm{x}}$ are statistically independent. With this approximation in mind, one can translate (3.35) into

Repeating the above argument then immediately yields

The case with $\lambda>0$ can be intuitively understood in a very similar way by first de-shrinking the estimate; we omit it here for brevity. We note that these hand-waving arguments can all be made rigorous, which is the main content of the proof.

8 Inference based on spectral estimates?

One would naturally be curious about whether there are other estimation procedures that also enable reasonable statistical inference. While this is beyond the scope of the current paper, we take a moment to discuss one alternative: the spectral method, as pioneered by [KMO10a, KMO10b] in the matrix completion problem. In a nutshell, this approach consists in computing a rank-$r$ approximation to $\mathcal{P}_{\Omega}(\bm{M})/p$, which is precisely the spectral initialization widely used in a two-stage nonconvex algorithm (cf. Algorithm 1) [KMO10b, SL16, CW15, CCF18, MWCC17]. While inference has not been, as far as we know, the focus of prior work on spectral methods,We note that inference from spectral estimates has been investigated in other context beyond matrix completion (e.g. the model without missing data [Xia19, FFHL19]). the recent papers [AFWZ17, MWCC17] hinted at the possibility of characterizing the distribution of the spectral estimate. Take a simple symmetric rank-1 case for example (i.e. $\bm{M}^{\star}=\bm{x}^{\star}\bm{x}^{\star\top}$ with $\|\bm{x}^{\star}\|_{2}=1$): the leading eigenvector $\bm{u}^{\mathsf{spectral}}$ of $\mathcal{P}_{\Omega}(\bm{M})/p$ often admits the following approximation (up to a global sign)

Expanding $\mathcal{P}_{\Omega}(\bm{M})=p\bm{x}^{\star}\bm{x}^{\star\top}+\mathcal{P}_{\Omega}(\bm{x}^{\star}\bm{x}^{\star\top})-p\bm{x}^{\star}\bm{x}^{\star\top}+\mathcal{P}_{\Omega}(\bm{E})$, we arrive at

In words, two major factors dictate the uncertainty of the spectral estimate: (1) the additive Gaussian noise (cf. the 1st term on the right-hand side of (3.42)), and (2) random sub-sampling (in particular, the randomness incurred by employing the sub-sampled $\mathcal{P}_{\Omega}(\bm{x}^{\star}\bm{x}^{\star\top})/p$ to approximate the truth $\bm{x}^{\star}\bm{x}^{\star\top}$). Given that the random sub-sampling effect cannot be ignored at all, the spectral estimates often suffer from a much larger estimation error (and hence a higher degree of uncertainty) compared to either the convex or the nonconvex estimates. In truth, this random sub-sampling effect does not go away even when the noise vanishes. Consequently, uncertainty quantification based on the spectral estimates may not be the most desirable option.

Prior art

Low-rank matrix completion, or more broadly, low-rank matrix recovery, is a fundamental task that permeates through a wide spectrum of applications in science, engineering, and finance (e.g. [RS05, SY07, CC14, FSZZ18, CCG15, ZPL15, BN06, CC18b, FWZ19, KS11, CZ16, KX15, FLM13, DR17, CDDD19, DPVW14, SZ12, FS11]). A paper of this length is unable to review all papers motivating and contributing to this enormous subject; interested readers are referred to [DR16, CC18a] for extensive discussions of motivating applications as well as the exciting recent development.

Numerous algorithms have been proposed to solve this problem efficiently, with two paradigms being arguably the most widely used: convex relaxation and nonconvex optimization. We briefly review the literature contributing to these two paradigms.

Nonconvex optimization algorithms, as pioneered by [KMO10a, Sre04], become increasingly more popular for solving various low-rank factorization problems, due to their appealing computational complexities [JNS13, CLS15, CC17, TBS+16, SL16, ZL16, CCFM19, WZG16, CLL19]. For instance, the gradient-based nonconvex methods have been analyzed for noisy matrix completion [KMO10b, CW15, MWCC17, CCF+19], which are shown to achieve near-optimal statistical accuracy and linear-time convergence guarantees all at once. Going beyond gradient methods, we note that other nonconvex methods (e.g. [RS05, JMD10, WYZ12, JNS13, FRW11, Van13, LXY13, Har14, JKN16, RT11, WCCL16, DC18, ZWL15, ZWYG18, MSL19, CCD+19]) and landscape properties [GLM16, CL17, GJZ17, ZSL19, ZJSL18, SXZ19] have been largely explored as well. The interested readers are referred to [CLC19] for an in-depth discussion. One limitation, however, is that the theoretical guarantees provided for nonconvex algorithms often exhibit sub-optimal dependency in the rank $r$ of the unknown matrix; for instance, most theory requires a sample complexity of at least $nr^{2}$ (in fact, often much larger than $nr^{2}$). This is outperformed by the convex relaxation approach.

Despite these recent developments, very little work has investigated statistical inference for noisy matrix completion. While [CKLN18, CKL16, CN15, CEGN15] discussed the construction of “honest” confidence regions, the volume of these regions is dependent on some (possibly huge) hidden constants, thus resulting in over-coverage. Perhaps the closest to our paper is the recent work [Xia18], which investigated inference for low-rank trace regression. Employing a closely related de-biased estimator with sample splitting, the paper [Xia18] established asymptotic normality of a certain projected distance between the estimate and the truth. The result therein, however, requires a sampling mechanism obeying the restricted isometry property (e.g. i.i.d. Gaussian designs), which fails to hold for matrix completion. Also, our approach does not require sample splitting — a technique that is convenient for analysis but conservative in constructing confidence regions. Another work by Cai et al. [CLR16] developed a unified approach to provide inference guarantees for linear inverse problems including low-rank matrix estimation. Their results, however, require the sample size to exceed the total dimension $n^{2}$ even under the Gaussian design. Finally, a recent line of work [MX17] explored uncertainty quantification under the Bayesian setting, hypothesizing on a special prior regarding the true matrix. This departs drastically from the scenario considered herein.

Inference in high-dimensional problems.

Inference in high-dimensional sparse regression has received much attention in the last few years [WR09, ZZ14, BCH11, vdGBRD14, JM14b, DBMM15, CG17, NL17, NNLL18, LSST16, LTTT14, MMB09, DBZ17, ZC17, BFL+18]. Our inferential approach is partly inspired by the recent developments on this topic, particularly with regard to the de-biased / de-sparsified estimators proposed for Lasso. More specifically, recognizing the non-negligible bias of the Lasso estimate

A line of work [ZZ14, vdGBRD14, JM14a] came up with a linear transformation of $\widehat{\bm{\beta}}$ of the form

Interestingly, our de-biased estimator (3.7) for matrix completion admits a very similar form as (4.2). To see this, recall that our de-biased estimator is given by

where $\bm{Z}$ can be either $\bm{Z}^{\mathsf{cvx}}$ or $\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}$ (see Table 1). Let $T$ be the tangent space of the set of rank-$r$ matrices at $\bm{Z}^{\mathsf{cvx},r}$ (resp. $\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}$) in the convex (resp. nonconvex) case, and $\mathcal{P}_{T}$ be the projection operator onto $T$. Somewhat surprisingly, replacing $\mathcal{P}_{\text{rank-}r}$ by $\mathcal{P}_{T}$ does not affect the de-biased estimator by much, in the sense that

In addition, recognizing that $\bm{Z}$ almost lies within the tangent space $T$,More precisely, if $\bm{Z}=\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}$, then $\bm{Z}\in T$; if $\bm{Z}=\bm{Z}^{\mathsf{cvx}}$, one has $\mathcal{P}_{T}(\bm{Z})\approx\bm{Z}$.one can rewrite

Finally, de-biased estimators have been put forward to tackle other high-dimensional problems, including but not limited to generalized linear models [vdGBRD14, NL17], graphical models [JVDG15, RSZZ15, MLL17, JvdG17], sparse PCA [JvdG18], treatment effects estimation [CCD+18, AIW18]. These are beyond the scope of the current paper.

Architecture of the proof

This section outlines the main steps for establishing Theorem 1 and Theorem 2. Before starting, we introduce some useful notation. For convenience of presentation, we insert the factor $1/p$ into (3.4) and redefine the nonconvex loss function as

In addition, for each $1\leq j,k\leq n$, we define the indicator $\delta_{jk}\triangleq\operatorname{\mathds{1}}\{(j,k)\in\Omega\}$, which is a Bernoulli random variable with mean $p$.

We also note that Theorem 1 (resp. Theorem 2) is subsumed by Theorem 5 (resp. Theorem 6). As a result, we shall focus on establishing Theorem 5 (resp. Theorem 6) when it comes to estimating low-rank factors (resp. the entries of the matrix).

Instate the assumptions of Theorem 5. Recall the definition of $v_{ij}^{\star}$ in (3.17). Then one has the following decomposition

where $g_{ij}\sim\mathcal{N}(0,v_{ij}^{\star})$ and the residual obeys — with probability exceeding $1-O(n^{-10})$ — that

Note that Theorem 5 and Theorem 6 are concerned with the de-biased estimators built upon both convex and nonconvex estimates. At first glance, one needs to establish theoretical guarantees for each of them separately. Fortunately, as alluded to previously (cf. (3.5)), the convex and nonconvex estimates are extremely close — a fact that has been established in [CCF+19]. The proximity of these two estimates naturally extends to the de-biased estimators constructed based on them. As a result, it suffices to concentrate on proving the theorems for any of these estimators; the claims for the other one follow immediately.

The following key lemma formalizes this argument, which will be established in Appendix C (see also Figure 4 for numerical evidence). Before continuing, we remind the readers of the key notation (see Appendix B for precise definitions):

$(\bm{X}^{\mathsf{ncvx}},\bm{Y}^{\mathsf{ncvx}})$: an approximate solution to the nonconvex problem (3.4) (see Appendix A.1);

$\bm{M}^{\mathsf{cvx,d}},\bm{X}^{\mathsf{cvx,d}},\bm{Y}^{\mathsf{cvx,d}}$: the de-biased estimators built upon the convex optimizer $\bm{Z}^{\mathsf{cvx}}$;

$\bm{M}^{\mathsf{ncvx,d}},\bm{X}^{\mathsf{ncvx,d}},\bm{Y}^{\mathsf{ncvx,d}}$: the de-biased estimators built upon the nonconvex estimate $(\bm{X}^{\mathsf{ncvx}},\bm{Y}^{\mathsf{ncvx}})$.

Suppose that the sample size obeys $n^{2}p\geq C\kappa^{4}\mu^{2}r^{2}n\log^{3}n$ for some sufficiently large constant $C>0$ and the noise satisfies $\sigma\sqrt{(\kappa^{4}\mu nr\log n)/p}\leq c\sigma_{\min}$ for some sufficiently small constant $c>0$. Set $\lambda=C_{\lambda}\sigma\sqrt{np}$ with some large enough constant $C_{\lambda}>0$.

With probability at least $1-O(n^{-10})$, one has

where $\mathcal{O}^{r\times r}$ is the set of $r\times r$ rotation matrices.

With probability exceeding $1-O(n^{-10})$, one has

In short, the first part of Lemma 3 tells us that

whereas the second part of Lemma 3 justifies that the proposed de-biased estimator is closely approximated by a linearized version (cf. (4.4)). Note that this linearized form bears a resemblance to the de-biased estimators developed for sparse linear regression [ZZ14, vdGBRD14, JM14a].

With Lemma 3 in place, we shall, from now on, focus on proving the main theorems for the nonconvex estimators, viz.

establishing Theorem 5 for the de-shrunken low-rank factors $(\bm{X}^{\mathsf{ncvx,d}},\bm{Y}^{\mathsf{ncvx,d}})$;

establishing Theorem 6 for the de-biased matrix estimator $\bm{X}^{\mathsf{ncvx,d}}\bm{Y}^{\mathsf{ncvx,d}\top}$.

To simplify the presentation hereafter, we shall use the following notation throughout the rest of this section:

$(\bm{X},\bm{Y})$: the nonconvex estimate $(\bm{X}^{\mathsf{ncvx}},\bm{Y}^{\mathsf{ncvx}})$;

$(\bm{X}^{\mathsf{d}},\bm{Y}^{\mathsf{d}})$: the de-shrunken estimate defined in (3.8) based on $(\bm{X},\bm{Y})=(\bm{X}^{\mathsf{ncvx}},\bm{Y}^{\mathsf{ncvx}})$;

$\bm{M}^{\mathsf{d}}\triangleq\bm{X}^{\mathsf{d}}\bm{Y}^{\mathsf{d}\top}$.

2 A precise characterization of the de-shrunken low-rank factors

We start with a precise characterization of the de-shrunken low-rank factors $\bm{X}^{\mathsf{d}}$ and $\bm{Y}^{\mathsf{d}}$, which paves the way for demonstrating both Theorem 5 and Theorem 6.

One has the following decompositions for $\bm{X}^{\mathsf{d}}$ and $\bm{Y}^{\mathsf{d}}$

which measures the imbalance between the low-rank factors $\bm{X}$ and $\bm{Y}$.

The claims follow from straightforward algebraic manipulations; see Appendix D.1. ∎

We make a few observations regarding Lemma 4. Take the decomposition of $\bm{X}^{\mathsf{d}}$ (5.8a) as an example:

First, the term $\bm{A}\bm{Y}^{\mathsf{d}}(\bm{Y}^{\mathsf{d}\top}\bm{Y}^{\mathsf{d}})^{-1}$ vanishes when we have full observations, i.e. $p=1$. Second, the terms involving $\nabla_{\bm{X}}f(\bm{X},\bm{Y})$ and $\bm{\Delta}_{\mathsf{balancing}}$ are both zero if $(\bm{X},\bm{Y})$ is an exact stationary point of $f(\cdot,\cdot)$; to see this, it is not hard to verify that any stationary point of $f(\cdot,\cdot)$ necessarily satisfies $\bm{X}^{\top}\bm{X}=\bm{Y}^{\top}\bm{Y}$, which in turn implies $\bm{\Delta}_{\mathsf{balancing}}=\bm{0}$. Consequently, the last three terms in (5.8a) are expected to be small when $p$ is sufficiently large and $(\bm{X},\bm{Y})$ is near a stationary point.

We shall make precise these arguments in subsequent subsections.

3 Taking global rotation into account

In order to invoke the decompositions of $\bm{X}^{\mathsf{d}}$ and $\bm{Y}^{\mathsf{d}}$ (cf. Lemma 4) to characterize the estimation errors, we still need to incorporate the (unrecoverable) rotation matrix. From now on, we shall focus primarily on the factor $\bm{X}^{\mathsf{d}}$. The claims on the other factor $\bm{Y}^{\mathsf{d}}$ can be easily obtained via symmetry.

where we recall that $\bm{H}^{\mathsf{d}}$ is the rotation matrix that best aligns $(\bm{X}^{\mathsf{d}},\bm{Y}^{\mathsf{d}})$ and $(\bm{X}^{\star},\bm{Y}^{\star})$ (see (3.12)). Substituting the identity

into the decomposition (5.8a), we arrive at

4 Key lemmas for establishing Theorem 5

We now state five key lemmas. Taking these collectively and substituting them into (5.11) immediately establish Theorem 5.

We shall start by controlling the term $\bm{\Phi}_{1}$ as defined in (5.12).

Suppose that the sample complexity obeys $n^{2}p\geq C\kappa^{4}\mu^{2}r^{2}n\log^{3}n$ for some sufficiently large constant $C>0$ and the noise satisfies $\sigma\sqrt{(\kappa^{4}\mu rn\log n)/p}\leq c\sigma_{\min}$ for some sufficiently small constant $c<0$. Then with probability at least $1-O(n^{-10})$, we have

Fix any $1\leq j\leq n$. If the de-shrunken estimate $\overline{\bm{Y}}^{\mathsf{d}}$ were independent of the randomness in the $j$th row of the matrix, i.e. $\bm{e}_{j}^{\top}\mathcal{P}_{\Omega}(\bm{E})$, then $\|\bm{e}_{j}^{\top}\bm{\Phi}_{1}\|_{2}$ would be well controlled. This hypothesis is certainly false, as $\overline{\bm{Y}}^{\mathsf{d}}$ clearly depends on $\bm{e}_{j}^{\top}\mathcal{P}_{\Omega}(\bm{E})$. Nevertheless, by exploiting the leave-one-out technique recently used in [EKBB+13, EK15, AFWZ17, MWCC17, CFMW19, CCF+19, CLL19, DC18], one can properly decouple the dependency and establish the desired bound. See Appendix D.2.∎

The next lemma controls the size of $\|\bm{\Phi}_{2}\|_{2,\infty}$. In essence, the term $\bm{\Phi}_{2}$ measures the difference between the estimate $\overline{\bm{Y}}^{\mathsf{d}}$ and the true signal $\bm{Y}^{\star}$; the closer these two are, the smaller $\|\bm{\Phi}_{2}\|_{2,\infty}$ should be. See Appendix D.3 for the proof of the following result.

Suppose that the sample complexity obeys $n^{2}p\geq C\kappa^{4}\mu^{2}r^{2}n\log^{3}n$ for some sufficiently large constant $C>0$ and the noise satisfies $\sigma\sqrt{(\kappa^{4}\mu rn\log n)/p}\leq c\sigma_{\min}$ for some sufficiently small constant $c<0$. Then with probability exceeding $1-O(n^{-10})$, one has

Moving on to $\bm{\Phi}_{3}$ and $\bm{\Phi}_{4}$, one has the following lemmas.

Suppose that the sample complexity obeys $n^{2}p\geq C\kappa^{4}\mu^{2}r^{2}n\log^{3}n$ for some sufficiently large constant $C>0$ and the noise satisfies $\sigma\sqrt{(\kappa^{4}\mu rn\log n)/p}\leq c\sigma_{\min}$ for some sufficiently small constant $c<0$. Then with probability exceeding $1-O(n^{-10})$, we have

It is straightforward to check that when $p=1$, one has $\|\bm{\Phi}_{3}\|_{2,\infty}=\|\bm{A}\|=0$, where we recall the definition of $\bm{A}$ in (5.7). Therefore, one expects $\|\bm{\Phi}_{3}\|_{2,\infty}$ to be small when $p$ is sufficiently large. See Appendix D.4.∎

Suppose that the sample complexity obeys $n^{2}p\geq C\kappa^{4}\mu^{2}r^{2}n\log^{3}n$ for some sufficiently large constant $C>0$ and the noise satisfies $\sigma\sqrt{(\kappa^{4}\mu rn\log n)/p}\leq c\sigma_{\min}$ for some sufficiently small constant $c<0$. Then with probability at least $1-O(n^{-10})$, one has

It is easily seen that the size of $\bm{\Phi}_{4}$ depends on how close $(\bm{X},\bm{Y})$ is to a stationary point of $f(\cdot,\cdot)$. For instance, in the extreme case where $(\bm{X},\bm{Y})$ is an exact stationary point, then one would have $\bm{\Phi}_{4}=\bm{0}$. See Appendix D.5.∎

The last lemma asserts that $\mathcal{P}_{\Omega}(\bm{E})\bm{Y}^{\star}(\bm{Y}^{\star\top}\bm{Y}^{\star})^{-1}/p$ is, in some sense, close to a zero-mean Gaussian random matrix with the desired covariance.

Suppose that the sample size obeys $n^{2}p\geq C\kappa^{2}\mu rn\log^{3}n$ for some sufficiently large constant $C>0$. Then one has the decomposition

Fix any $1\leq j\leq n$. The $j$th row, namely $\bm{e}_{j}^{\top}[\mathcal{P}_{\Omega}(\bm{E})\bm{Y}^{\star}(\bm{Y}^{\star\top}\bm{Y}^{\star})^{-1}/p]$ is conditionally Gaussian in the sense that

where we recall that $\delta_{jk}=\operatorname{\mathds{1}}\{(j,k)\in\Omega\}$. Recognize that the conditional covariance matrix concentrates sharply around its expectation, i.e. $\sigma^{2}(\bm{\Sigma}^{\star})^{-1}/p$, which is the covariance matrix of $\bm{Z}_{\bm{X}}^{\top}\bm{e}_{j}$ that we are after. Hence, one can expect that $\mathcal{P}_{\Omega}(\bm{E})\bm{Y}^{\star}(\bm{Y}^{\star\top}\bm{Y}^{\star})^{-1}/p$ is, marginally, not too far from a Gaussian random matrix. This argument can be carried out formally; see Appendix D.6. ∎

5 From low-rank factors to matrix entries (Proof of Theorem 6)

We now turn attention to inference on the matrix entries, by establishing Theorem 6. Towards this, we first make the following observation: for any $1\leq i,j\leq n$,

One can readily apply the decompositions in Theorem 5 to obtain

Take the preceding three identities collectively to reach

Following the same route as in Section 5.4, one can verify that $\Theta_{ij}=\bm{e}_{i}^{\top}\bm{Z}_{\bm{X}}\bm{Y}^{\star\top}\bm{e}_{j}+\bm{e}_{i}^{\top}\bm{X}^{\star}\bm{Z}_{\bm{Y}}^{\top}\bm{e}_{j}$ is approximately Gaussian, whereas the residual term $\Lambda_{ij}$ is small in magnitude. These claims are formally stated in the next two lemmas, with the proofs deferred to Appendix E.

Suppose that the sample complexity obeys $n^{2}p\geq C\kappa^{4}\mu^{2}r^{2}n\log^{3}n$ for some sufficiently large constant $C>0$ and the noise satisfies $\sigma\sqrt{(\kappa^{4}\mu rn\log n)/p}\leq c\sigma_{\min}$ for some sufficiently small constant $c<0$. Then with probability exceeding $1-O(n^{-10})$, one has

Suppose that $np\geq C\kappa^{2}\mu r^{2}\log^{2}n$ for some sufficiently large constant $C>0$. Then we have the decomposition

where $g_{ij}\sim\mathcal{N}(0,v_{ij}^{\star})$ and the remaining term $\theta_{ij}$ satisfies — with probability exceeding $1-O(n^{-10})$ — that

Putting the above two lemmas together immediately establishes Theorem 6 and hence Theorem 2.

Discussion

The present paper makes progress towards inference and uncertainty quantification for noisy matrix completion, by developing simple de-biased estimators that admit tractable and accurate distributional characterizations. While we have achieved some early success in accomplishing this, our results are likely sub-optimal in the following aspects:

Dependency on the rank and the condition number. To enable valid inference, our sample complexity (cf. (3.13) and (3.18a)) scales sub-optimally with the rank $r$ and the condition number $\kappa$. The sub-optimality can be understood through comparisons with the sample size requirement $O(nr\log^{2}n)$ in the noise-free settings, which is independent of $\kappa$ and matches the information limit (up to some $\log$ factor). Improving such dependency calls for more refined analysis techniques.

Detection of the size of the entries. On one hand, when the size of the entry $M_{ij}^{\star}$ is moderately large (cf. (3.18b)), Corollary 1 allows us to construct a valid confidence interval for it. On the other hand, when $\|\bm{U}_{i,\cdot}^{\star}\|_{2}+\|\bm{V}_{j,\cdot}^{\star}\|_{2}$ vanishes, Theorem 5 tells us that the estimation error $M_{ij}^{\mathsf{d}}-M_{ij}^{\star}$ is better approximated by the inner product of two independent Gaussian random vectors. It remains to be seen how to determine whether $\|\bm{U}_{i,\cdot}^{\star}\|_{2}+\|\bm{V}_{j,\cdot}^{\star}\|_{2}$ is too small.

Low signal-to-noise (SNR) regime. Our theory operates under the moderate-to-high SNR regime, where $\sigma^{2}_{\min}/\sigma^{2}$ (which is proportional to the SNR) is required to exceed the order of $n/p$; see the conditions in Theorem 5. It is unclear whether the connection between the convex and the nonconvex estimators hold in the low SNR regime. How to conduct inference in such a scenario is an important future direction.

In addition, our investigation has been dedicated to a natural random model, which by no means covers the most general settings of practical interest. There are numerous possible extensions that merit future investigation:

Approximate low-rank structure. Our current theory is built upon the exact low-rank structure of $\bm{{M}}^{\star}$. Realistically, the matrix of interest is often only approximately low-rank. It is of great interest to study how to carry out statistical inference under such imperfect structural assumptions.

More general sampling patterns. This paper operates under the uniform random sampling assumption, which might sometimes be off in practical situations. It would be interesting to investigate whether our results in this paper can extend to more general non-uniform sampling patterns (e.g. [NW12]).

Extensions to robust PCA, sparse PCA, and 1-bit matrix completion. A variety of important extensions of matrix completion have been explored in prior literature, including but not limited to the case with sparse outliers (i.e. robust PCA [CLMW11, CSPW11]), the case where the matrix of interest is simultaneously sparse and low-rank (i.e. sparse PCA [ZHT06, CMW13]), and the case where only finite-bit observations are available (i.e. 1-bit matrix completion [DPVW14, CZ13]). Performing valid uncertainty assessment for these scenarios requires non-trivial extensions of the link between convex and nonconvex optimization.

Other loss functions. In the estimation stage, one might sometimes prefer other loss functions beyond the penalized squared loss. This might arise from either statistical considerations (e.g. employing a penalized Poisson log-likelihood to accommodate Poisson noise [CX16]), or computational concerns (e.g. adopting a non-smooth loss to improve convergence [CCD+19]). It would be of fundamental importance to develop a unified inferential framework that covers a broader family of loss functions.

Acknowledgements

Y. Chen is supported in part by the AFOSR YIP award FA9550-19-1-0030, by the ONR grant N00014-19-1-2120, by the ARO grant W911NF-18-1-0303, by the NSF grants CCF-1907661 and IIS-1900140, and by the Princeton SEAS innovation award. J. Fan is supported in part by NSF grants DMS-1662139 and DMS-1712591, ONR grant N00014-19-1-2120, and NIH grant 2R01-GM072611-13. C. Ma is supported in part by Hudson River Trading AI Labs (HAIL) Fellowship. This work was done in part while Y. Chen was visiting the Kavli Institute for Theoretical Physics (supported in part by the NSF grant PHY-1748958). We thank Weijie Su for helpful discussions.

Appendix A Preliminaries

In this section, we gather several notation and preliminary facts that are useful throughout the analysis. All the proofs, if needed, are deferred to Appendix I.

To begin with, we make precise the algorithm used to minimize the nonconvex loss function (5.1). Specifically, we describe the following details that are crucial for us to implement Algorithm 1:

Set the initial point to be $(\bm{X}^{0},\bm{Y}^{0})=(\bm{X}^{\star},\bm{Y}^{\star})$ or the spectral initialization as in [MWCC17, CLL19];

Set the stepsize $\eta\asymp 1/(n^{6}\kappa^{3}\sigma_{\max})$;

Set the maximum number of iterations to be $t_{0}\asymp n^{23}$;

The returned estimate is $(\bm{X}^{\mathsf{ncvx}},\bm{Y}^{\mathsf{ncvx}})\triangleq(\bm{X}^{t_{\star}},\bm{Y}^{t_{\star}})$, where

In words, we run gradient descent in Algorithm 1 until we reach a point whose gradient is exceedingly small.

Many of the preliminary facts below were established for the case $(\bm{X}^{0},\bm{Y}^{0})=(\bm{X}^{\star},\bm{Y}^{\star})$ [CCF+19], which is certainly not implementable in practice, however, it serves as a good proxy for studying the convex estimator. Fortunately, the same theoretical guarantees stated in Appendix A.2 can be readily established for spectral initialization using almost the same arguments adopted in [MWCC17, CLL19, CCF+19]. We omit this part mainly for the sake of brevity.

To facilitate analysis, we introduce a set of auxiliary nonconvex loss functions. For any $1\leq j\leq n$, define

be the nonconvex estimate returned by this auxiliary algorithm (i.e. Algorithm 2), which serves as an approximate solution to (A.2).

A.2 Properties of approximate nonconvex solutions

This subsection gathers the properties of the (approximate) nonconvex solutions. Throughout this subsection, we use the shorthand

and recall the definition of $(\bm{X}^{(j)},\bm{Y}^{(j)})$ in (A.3). The regularization parameter is chosen to satisfy

The claims stated below hold under the sample complexity and the noise condition presumed in [CCF+19, Theorem 1] (see also Theorem 5 in the current manuscript)

The first set of facts is related to $(\bm{X},\bm{Y})$. In view of [CCF+19], $\bm{F}$ is a faithful estimateTechnically, the statements in [CCF+19, Lemma 5] are for $\eta\asymp 1/(n\kappa^{3}\sigma_{\max})$ and $t_{0}\asymp n^{18}$. Nevertheless, inspecting their proofs reveals that the claims continue to hold for our choices $\eta\asymp 1/(n^{6}\kappa^{3}\sigma_{\max})$ and $t_{0}\asymp n^{23}$. of $\bm{F}^{\star}$, in the sense that

hold with probability exceeding $1-O(n^{-10})$. In addition, on the same high-probability event, one has

In words, the first claim ensures that $(\bm{X},\bm{Y})$ is an approximate stationary point of $f(\cdot,\cdot)$; the second bound tells us that $(\bm{X},\bm{Y})$ is nearly balanced, in the sense that $\bm{X}^{\top}\bm{X}\approx\bm{Y}^{\top}\bm{Y}$; the last one formalizes the proximity between the convex solution and the nonconvex one; see also (3.5).

hold with probability at least $1-O(n^{-10})$.

As has been shown in [CCF+19], the leave-one-out auxiliary point $(\bm{X}^{(j)},\bm{Y}^{(j)})$ satisfies

with probability exceeding $1-O(n^{-10})$.

Parallel to the transition from $(\bm{X},\bm{Y})$ to $(\bm{X}^{\mathsf{d}},\bm{Y}^{\mathsf{d}})$, we set

to be the de-shrunken estimators of $\bm{X}^{(j)}$ and $\bm{Y}^{(j)}$, respectively. We shall demonstrate in Appendix I that, with probability at least $1-O(n^{-10})$,

In addition to these four sets of claims, we have the following immediate consequence of the incoherence condition (2.4)

Moreover, recall that $\bm{A}=(1/p)\cdot\mathcal{P}_{\Omega}\left(\bm{X}\bm{Y}^{\top}-\bm{X}^{\star}\bm{Y}^{\star\top}\right)-\left(\bm{X}\bm{Y}^{\top}-\bm{X}^{\star}\bm{Y}^{\star\top}\right)$ (cf. (5.7)). We obtain from the proof of [CCF+19, Lemma 8] that

Last but not least, we list a few simple but useful results: the nonconvex solution $\bm{F}$ satisfies

The same holds true if we replace $\bm{F}$ by either $\bm{F}^{\mathsf{d}}$, $\bm{F}^{(j)}$ $\bm{F}^{\mathsf{d},(j)}$ or their corresponding low-rank factors. Here $j$ can vary from 1 to $n$.

Appendix B Summary of the proposed estimators

Let $\bm{Z}^{\mathsf{cvx}}$ be the minimizer of the convex program (3.1), and let $(\bm{X}^{\mathsf{ncvx}},\bm{Y}^{\mathsf{ncvx}})$ be the solution returned by the Algorithm 1 (with algorithmic details specified in Appendix A.1). Recall that $\bm{Z}^{\mathsf{cvx},r}$ is the best rank-$r$ approximation of $\bm{Z}^{\mathsf{cvx}}$, viz.

In addition, we let the matrix estimate obtained by the nonconvex algorithm be $\bm{Z}^{\mathsf{ncvx}}\triangleq\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}$. We further denote by $(\bm{X}^{\mathsf{cvx}},\bm{Y}^{\mathsf{cvx}})$ the estimate of low-rank factors obtained by convex relaxation; more specifically, we set $(\bm{X}^{\mathsf{cvx}},\bm{Y}^{\mathsf{cvx}})$ to be the balanced rank-$r$ factorization of $\bm{Z}^{\mathsf{cvx},r}$ obeying $\bm{X}^{\mathsf{cvx}}\bm{Y}^{\mathsf{cvx}\top}=\bm{Z}^{\mathsf{cvx},r}$ and $\bm{X}^{\mathsf{cvx}\top}\bm{X}^{\mathsf{cvx}}=\bm{Y}^{\mathsf{cvx}\top}\bm{Y}^{\mathsf{cvx}}$. With these notations in place, our de-biased and de-shrunken estimators can be summarized as follows.

De-shrunken estimators for low-rank factors:

Appendix C Proof of Lemma 3

Throughout this section, let $\bm{U}\bm{\Sigma}\bm{V}^{\top}$ be the rank-$r$ SVD of the nonconvex estimate $\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}$ and $T$ the tangent space of the set of rank-$r$ matrices at $\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}$. Correspondingly, we denote by $\mathcal{P}_{T}$ the projection operator onto the tangent space $T$, and let $\mathcal{P}_{T^{\perp}}=\mathcal{I}-\mathcal{P}_{T}$, where $\mathcal{I}$ is the identity operator.

In essence, we intend to justify that $\bm{M}^{\mathsf{cvx,d}}$, $\bm{M}^{\mathsf{ncvx,d}}$ and $\bm{X}^{\mathsf{ncvx,d}}\bm{Y}^{\mathsf{ncvx,d}\top}$ are all very close to $\bm{U}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})\bm{V}^{\top}$.

Recall from the definition of the de-biased estimator $\bm{M}^{\mathsf{cvx,d}}$ (cf. (B.1a)) that

Replacing $\bm{Z}^{\mathsf{cvx}}$ by $\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}$ results in

Apply the proximity bound (A.12) to obtain (recall that in (A.12), one has $(\bm{X},\bm{Y})=(\bm{X}^{\mathsf{ncvx}},\bm{Y}^{\mathsf{ncvx}})$)

as long as $n^{5}p\gg\kappa^{2}$. In addition, in view of [CCF+19, Claim 2], one has the decomposition

where the middle line follows since $\bm{U}\bm{\Sigma}\bm{V}^{\top}$ is defined to be the SVD of $\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}$.

where the last inequality results from (C.3) and (C.5). Combining the above two bounds with the fact that $\bm{U}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})\bm{V}^{\top}$ and $\mathcal{P}_{T}(\bm{\Delta}_{\bm{Z}}-\frac{1}{p}\bm{R})$ are orthogonal to each other, we arrive at the conclusion that $\bm{U}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})\bm{V}^{\top}$ is the top-$r$ SVD of $\bm{C}$ and

Second, let $\hat{\bm{U}}\hat{\bm{\Sigma}}\hat{\bm{V}}^{\top}$ be the top-$r$ SVD of $\bm{C}+\bm{\Delta}$. By definition, one has $\hat{\bm{U}}\hat{\bm{\Sigma}}\hat{\bm{V}}^{\top}=\bm{M}^{\mathsf{cvx,d}}$. We are left with checking the two conditions in Lemma 14. To begin with, the perturbation term $\bm{\Delta}$ obeys

where (i) comes from (C.3) and (C.5) and the last inequality (ii) arises since $np\gg\kappa^{2}.$ Clearly, the size of the perturbation is much smaller than $\lambda/p$ and hence $\|\bm{C}\|$ (cf. (C.8a)). In addition,

where the equality depends on (C.8b) and the last line results from (C.7) and (C.9). Consequently,

Here the first relation arises from (C.8a) and the second holds since $\sigma_{r}(\bm{\Sigma})\geq\sigma_{\min}/2$, a simple consequence of (A.19). We are now ready to apply Lemma 14 to obtain

where we have used the fact that $\|\bm{\Sigma}+(\lambda/p)\bm{I}_{r}\|\lesssim\sigma_{\max}$, which also can be derived from (A.19). The above bound combined with (C.9) yields

as long as $np\gg\kappa^{3}$. We remark that by setting $\bm{\Delta}_{\bm{Z}}=\bm{0}$, one also obtains the bound on $\bm{M}^{\mathsf{ncvx,d}}$, i.e.

We move on to investigating $\|\bm{X}^{\mathsf{ncvx,d}}\bm{Y}^{\mathsf{ncvx,d}\top}-\bm{U}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})\bm{V}^{\top}\|$, for which we have the following claim.

Taking the above three bounds collectively and recognizing that $\lambda\lesssim\sigma\sqrt{np}$ yield the advertised bound (5.4a).

The last inequality is the small-gradient condition (see (A.10), in which $(\bm{X},\bm{Y})=(\bm{X}^{\mathsf{ncvx}},\bm{Y}^{\mathsf{ncvx}}))$. Employ the definitions for $\bm{X}^{\mathsf{ncvx,d}}$ and $\bm{Y}^{\mathsf{ncvx,d}}$ (cf. (B.2a) and (B.2b)) to see that

It then boils down to showing that (i) $\bm{{A}}_{1}$ is very close to $\bm{U}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})\bm{V}^{\top}$, and (ii) $\bm{{A}}_{2}$ is small in size.

First, recall that $\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}=\bm{{U}}\bm{{\Sigma}}\bm{{V}}^{\top}$, which combined with (C.12) gives

Here, the last inequality comes from (C.13) and its immediate consequence that $\|\bm{{\Sigma}}_{\bm{{Q}}}\|\asymp\|\bm{{\Sigma}}_{\bm{{Q}}}^{-1}\|\asymp 1$.

Second, apply the perturbation bound for matrix square roots (cf. Lemma 13) to obtain

Here, the inequality (i) depends on the facts that

whereas the inequality (ii) holds because of the facts that $\|(\bm{X}^{\mathsf{ncvx}\top}\bm{X}^{\mathsf{ncvx}})^{-1}\|\lesssim 1/\sigma_{\min}$, $\|(\bm{Y}^{\mathsf{ncvx}\top}\bm{Y}^{\mathsf{ncvx}})^{-1}\|\lesssim 1/\sigma_{\min}$ and the balancedness condition (A.11)

As a result, the operator norm of $\bm{{A}}_{2}$ is bounded by

Take (C.14), (C.15) and (C.17) collectively to arrive at

C.2 Proof of the inequality (5.4b)

Next, we switch attention to the low-rank factors. Our goal is to demonstrate that $(\bm{X}^{\mathsf{cvx,d}},\bm{Y}^{\mathsf{cvx,d}})$ and $(\bm{X}^{\mathsf{ncvx,d}},\bm{Y}^{\mathsf{ncvx,d}})$ are both extremely close to $(\bm{U}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})^{1/2},\bm{V}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})^{1/2})$ modulo some global rotation, which will be established in (C.19) and (C.20) shortly.

We start by justifying the proximity between $(\bm{X}^{\mathsf{ncvx,d}},\bm{Y}^{\mathsf{ncvx,d}})$ and $(\bm{U}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})^{1/2},\bm{V}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})^{1/2})$. In view of (C.12), we know that

The latter bound follows from similar derivations as in (C.16). A similar bound holds for $\bm{Y}^{\mathsf{ncvx,d}}$. Recognizing that

Next, we establish the connection between $(\bm{X}^{\mathsf{cvx,d}},\bm{Y}^{\mathsf{cvx,d}})$ and $(\bm{U}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})^{1/2},\bm{V}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})^{1/2})$. To accomplish this, we first study the relationship between $(\bm{X}^{\mathsf{cvx}},\bm{Y}^{\mathsf{cvx}})$ and $(\bm{U}\bm{\Sigma}^{1/2},\bm{V}\bm{\Sigma}^{1/2})$. Recall that $\bm{X}^{\mathsf{cvx}}$ and $\bm{Y}^{\mathsf{cvx}}$ constitute a balanced factorization of $\bm{Z}^{\mathsf{cvx},r}$, while $(\bm{U}\bm{\Sigma}^{1/2},\bm{V}\bm{\Sigma}^{1/2})$ is a balanced one of $\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}=\bm{U}\bm{\Sigma}\bm{V}^{\top}$. Hence one can view $\bm{Z}^{\mathsf{cvx},r}$ as a perturbation of $\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}=\bm{U}\bm{\Sigma}\bm{V}^{\top}$ and investigate the perturbation bounds on the balanced factorizations. Going through the same derivations as in [MWCC17, Appendix B.7] and [CLL19, Appendix B.2.1], one reaches

Here the last relation follows from the proximity of the convex estimator and the nonconvex estimator; see (A.12). Repeating the same argument as above to translate the bound between $(\bm{X}^{\mathsf{cvx}},\bm{Y}^{\mathsf{cvx}})$ and $(\bm{U}\bm{\Sigma}^{1/2},\bm{V}\bm{\Sigma}^{1/2})$ to that of $(\bm{X}^{\mathsf{cvx,d}},\bm{Y}^{\mathsf{cvx,d}})$ and $(\bm{U}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})^{1/2},\bm{V}(\bm{\Sigma}+\frac{\lambda}{p}\bm{I}_{r})^{1/2})$, we conclude that

This together with (C.19) and the assumption $n\gg\kappa^{4}$ concludes the proof.

C.3 Proof of the inequality (5.5)

We shall focus on proving the claim for the nonconvex estimator $\bm{M}^{\mathsf{ncvx,d}}$ and $\bm{X}^{\mathsf{ncvx}}\bm{Y}^{\mathsf{ncvx}\top}$; the claim for the convex estimator $\bm{Z}^{\mathsf{cvx}}$ can be treated similarly.