Bridging Convex and Nonconvex Optimization in Robust PCA: Noise, Outliers, and Missing Data

Introduction

A diverse array of science and engineering applications (e.g. video surveillance, joint shape matching, graph clustering, covariance modeling, graphical models) involves estimation of low-rank matrices [CLC19, CLMW11, CGH14, JCSX11, CPW12, FLM13, DR16]. The imperfectness of data acquisition processes, however, presents several common yet critical challenges: (1) random noise: which reflects the uncertainty of the environment and/or the measurement processes; (2) outliers: which represent a sort of corruption that exhibits abnormal behavior; and (3) incomplete data, namely, we might only get to observe a fraction of the matrix entries. Low-rank matrix estimation algorithms aimed at addressing these challenges have been extensively studied under the umbrella of robust principal component analysis (Robust PCA) [CSPW11, CLMW11], a terminology popularized by the seminal work [CLMW11].

where $\bm{S}^{\star}=[S_{ij}^{\star}]$ is a matrix consisting of outliers, $\bm{E}=[E_{ij}]$ represents the random noise, and we only observe entries over an index subset $\Omega_{\mathsf{obs}}\subseteq[n_{1}]\times[n_{2}]$ with $[n]\coloneqq\{1,2,\cdots,n\}$. The current paper assumes that $\bm{S}^{\star}$ is a relatively sparse matrix whose non-zero entries might have arbitrary magnitudes. This assumption has been commonly adopted in prior work to model gross outliers, while enabling reliable disentanglement of the outlier component and the low-rank component [CSPW11, CLMW11, CJSC13, Li13]. In addition, we suppose that the entries $\{E_{ij}\}$ are independent zero-mean sub-Gaussian random variables, as commonly assumed in the statistics literature to model a large family of random noise. The aim is to reliably estimate $\bm{L}^{\star}$ given the grossly corrupted and possibly incomplete data (1.1). Ideally, this task should be accomplished without knowing the locations and magnitudes of the outliers $\bm{S}^{\star}$.

Focusing on the noiseless case with $\bm{E}=\bm{0}$, the papers by [CSPW11, CLMW11] delivered a positive and somewhat surprising message: both the low-rank component $\bm{L}^{\star}$ and the sparse component $\bm{S}^{\star}$ can be efficiently recovered with absolutely no error by means of a principled convex program

Moving on to the more realistic noisy setting, a natural strategy is to solve the following regularized least-squares problem

Where does the algorithm (1.3) stand in terms of its statistical performance vis-à-vis random noise, sparse outliers and missing data?

Unfortunately, however simple this program (1.3) might seem, the existing theoretical support remains far from satisfactory, as we shall discuss momentarily.

2 Theory-practice gaps under random noise

To assess the tightness of prior statistical guarantees for (1.3), we find it convenient to first look at a simple setting where (i) $n_{1}=n_{2}=n$, (ii) $\bm{E}$ consists of independent Gaussian components, namely, $E_{ij}\sim\mathcal{N}(0,\sigma^{2})$, and (iii) there is no missing data. This simple scenario is sufficient to illustrate the sub-optimality of prior theory.

The paper [ZLW+10] was the first to derive a sort of statistical performance guarantees for the above convex program. Under mild conditions, [ZLW+10] demonstrated that any minimizer $(\widehat{\bm{L}},\widehat{\bm{S}})$ of (1.3) achievesMathematically, [ZLW+10] investigated an equivalent constrained form of (1.3) and developed an upper bound on the corresponding estimation error.

Consider an idealistic scenario where an oracle informs us of the outlier matrix $\bm{S}^{\star}$. With the assistance of this oracle, the task of estimating $\bm{L}^{\star}$ reduces to a low-rank matrix denoising problem [DG14]. By fixing $\bm{S}$ to be $\bm{S}^{\star}$ in (1.3), we arrive at a simplified convex program

It is known that (e.g. [DG14, CCF+20]): under mild conditions and with a properly chosen $\lambda$, the estimation error of (1.5) satisfies

where we abuse the notation and denote by $\widehat{\bm{L}}$ the minimizer of (1.5). The large gap between the above two bounds (1.4) and (1.6) is self-evident; in particular, if $r=O(1)$, the gap between these two bounds can be as large as an order of $n^{1.5}$.

All in all, there seems to be a large gap between the practical performance of (1.3) and the existing theoretical support. This calls for a new theory that better explains practice, which we pursue in the current paper. We remark in passing that statistical guarantees have been developed in [ANW12, KLT17] for other convex estimators (i.e. the ones that are different from the convex estimator (1.3) considered herein). We shall compare our results with theirs later in Section 1.4.

3 Models, assumptions and notation

represent the singular values and the condition number of $\bm{L}^{\star}$, respectively. We denote by $\Omega^{\star}$ the support set of $\bm{S}^{\star}$, that is,

With this set of notation in place, we list below our key model assumptions.

The low-rank matrix $\bm{L}^{\star}$ with SVD $\bm{L}^{\star}=\bm{U}^{\star}\bm{\Sigma}^{\star}\bm{V}^{\star\top}$ is assumed to be $\mu$-incoherent in the sense that

Each entry is observed independently with probability $p$, namely,

Each observed entry is independently corrupted by an outlier with probability $\rho_{\mathsf{s}}$, namely,

where $\Omega^{\star}\subseteq\Omega_{\mathsf{obs}}$ is the support of the outlier matrix $\bm{S}^{\star}$.

The signs of the nonzero entries of $\bm{S}^{\star}$ are i.i.d. symmetric Bernoulli random variables (independent from the locations), namely,

The noise matrix $\bm{E}=[E_{ij}]$ is composed of independent symmetricIn fact, we only require $E_{ij}$ to be symmetric for all $(i,j)\in\Omega^{\star}$. zero-mean sub-Gaussian random variables with sub-Gaussian norm at most $\sigma>0$, i.e. $\|E_{ij}\|_{\psi_{2}}\leq\sigma$ (see [Ver12, Definition 5.7] for precise definitions).

We take a moment to expand on our model assumptions. Assumption 1 is standard in the low-rank matrix recovery literature [CR09, CLMW11, Che15, CLC19]. If $\mu$ is small, then this assumption specifies that the singular spaces of $\bm{L}^{\star}$ is not sparse in the standard basis, thus ensuring that $\bm{L}^{\star}$ is not simultaneously low-rank and sparse. Assumption 3 requires the sparsity pattern of the outliers $\bm{S}^{\star}$ to be random, which precludes it from being simultaneously sparse and low-rank. In essence, Assumptions 1 and 3 are identifiability conditions, taken together as a sort of separation condition on $(\bm{L}^{\star},\bm{S}^{\star})$, which plays a crucial role in guaranteeing exact recovery in the noiseless case (i.e. $\bm{E}=\bm{0}$); see [CLMW11] for more discussions on these conditions. Assumption 4 requires the signs of the outliers to be random, which has also been made in [ZLW+10, WL17].Note that while the theorems in [ZLW+10, WL17] do not make explicit this random sign assumption, the proofs therein do rely on this assumption to guarantee the existence of certain approximate dual certificates. We shall discuss in detail the crucial role of this random sign assumption (as opposed to deterministic sign patterns) in Section 1.6.

4 Main results

Armed with the above model assumptions, we are positioned to present our improved statistical guarantees for convex relaxation (1.3) in the random noise setting. Without loss of generality, assume that

As we shall elucidate in Section 1.5 and Section 3, our theory is established by exploiting an intriguing and intimate connection between convex relaxation and nonconvex optimization, and hence the title of this paper.

For the sake of simplicity, we shall start by presenting our statistical guarantees when the rank $r$, the condition number $\kappa$ and the incoherence parameter $\mu$ of $\bm{L}^{\star}$ are all bounded by some constants. Despite its simplicity, this setting subsumes as special cases a wide array of fundamentally important applications, including angular and phase synchronization [Sin11] in computational biology, joint shape mapping problem [HG13, CGH14] in computer vision, and so on. All of these problems involve estimating a very well-conditioned matrix $\bm{L}^{\star}$ with a small rank.

Suppose that Assumptions 1-5 hold, and that $r,\kappa,\mu=O(1)$. Take $\lambda=C_{\lambda}\sigma\sqrt{n_{1}p}$ and $\tau=C_{\tau}\sigma\sqrt{\log n_{2}}$ in (1.3) for some large enough constants $C_{\lambda},C_{\tau}>0$. Assume that

for some sufficiently large constant $C_{\mathsf{sample}}>0$ and some sufficiently small constants $c_{\mathsf{noise}},c_{\mathsf{outlier}}>0$. Then with probability exceeding $1-O(n_{2}^{-3})$, the following holds:

Any minimizer $(\bm{L}_{\mathsf{cvx}},\bm{S}_{\mathsf{cvx}})$ of the convex program (1.3) obeys

and the statistical guarantees (1.13) hold unchanged if $\bm{L}_{\mathsf{cvx}}$ is replaced by $\bm{L}_{\mathsf{cvx},r}$.

Before we embark on interpreting our statistical guarantees, let us first parse the required conditions (1.12) in Theorem 1. For simplicity we assume that $n_{1}=n_{2}=n$.

Noise levels. The noise condition, namely $\sigma\sqrt{n\log n/p}\lesssim\sigma_{\min}$, accommodates a wide range of noise levels. To see this, it is straightforward to check that this noise condition is equivalent to

as long as $r,\mu,\kappa\asymp 1$. In other words, the entrywise noise level $\sigma$ is allowed to be significantly larger than the maximum magnitude of the entries in the low-rank matrix $\bm{L}^{\star}$, as long as $p\gg(\log n)/n$.

Tolerable fraction of outliers. The above theorem assumes that no more than a fraction $\rho_{s}\lesssim 1/\log n$ of observations are corrupted by outliers. In words, our theory allows nearly a constant proportion (up to a logarithmic order) of the entries of $\bm{L}^{\star}$ to be corrupted with arbitrary magnitudes.

Next, we move on to the interpretation of our statistical guarantees. Note that we still assume that $n_{1}=n_{2}=n$ for ease of presentation.

Near-optimal statistical guarantees. Our first result (1.13a) gives an Euclidean estimation error bound of (1.3)

with the proviso that $p$ is at least on the constant order. The restriction on $p$ arises from the dual certificate constructed in [CLMW11], which is also used in the Proof of Theorem 4 in [WL17]. While this is sub-optimal compared to our results in the setting considered herein, it is worth pointing out that the bound therein accommodates arbitrary noise matrix $\bm{E}$ (e.g. deterministic, adversary), and here in (1.16) we specialize their result to the random noise setting, namely the noise $\bm{E}$ obeys Assumption 5. In addition, under the full observation (i.e. $p=1$) setting, the paper [ANW12] derived an estimation error bound for a convex program similar to (1.3), but with an additional constraint regularizing the spikiness of the low-rank component. Note that instead of imposing the incoherence condition as in Assumption 1, the prior work [ANW12] assumes a milder spikiness condition on $\bm{L}^{\star}$, which only constrains the maximum entry in the matrix $\bm{L}^{\star}$ is not too large. When $\{E_{ij}\}$ are i.i.d. drawn from $\mathcal{N}(0,\sigma^{2})$ and when there is no missing data (i.e. $p=1$), the Euclidean estimation error bound achievable by their estimator $\bm{L}_{\mathsf{cvx}}^{\mathsf{ANW}}$ reads

which is sub-optimal compared to our results. In particular, (i) the bound (1.17) does not vanish even as the noise level decreases to zero, and (ii) it becomes looser as $\rho_{\mathsf{s}}$ grows (e.g. if $\rho_{\mathsf{s}}\asymp 1/\log n$, the bound (1.17) is $O(\sqrt{n})$ larger than our bound). Moreover, the work [ANW12] did not account for missing data. Similar to [ANW12] (but with an additional spikiness condition on $\bm{S}^{\star}$), the paper [KLT17] derived an estimation error bound for a constrained convex program, with a new constraint regularizing the spikiness of the sparse outliers. Their Euclidean estimation error bound reads

which is also sub-optimal compared to our results. In particular, (1) their error bound degrades as the magnitude $\|\bm{S}^{\star}\|_{\infty}$ of the outlier increases; (2) when there is no missing data (i.e. $p=1$), their bound might be off by a factor as large as $O(\sqrt{n})$. It is worth emphasizing that the theory developed in these prior works is developed to accommodate a broader range of matrices. For example, both [ANW12] and [KLT17] study the set of entrywise bounded low-rank matrices (without assuming the incoherence condition); [ANW12] even allows $\bm{L}^{\star}$ to be approximately low rank. To ease comparison, Table 1 displays a summary of our results vs. prior statistical guarantees when specialized to the settings considered herein.

as long as $r,\kappa,\mu\asymp 1$, which is about $O(n)$ times small than the Euclidean loss (1.15) modulo some logarithmic factor. This uncovers an appealing “delocalization” behavior of the estimation errors, namely, the estimation errors of $\bm{L}^{\star}$ are fairly spread out across all entries. This can also be viewed as an “implicit regularization” phenomenon: the convex program automatically controls the spikiness of the low-rank solution, without the need of explicitly regularizing it (e.g. adding a constraint $\|\bm{L}\|_{\infty}\leq\alpha$ as adopted in the prior work [ANW12, KLT17]). See Figure 2 for the numerical evidence for the relative entrywise and spectral norm error of $\bm{L}_{\mathsf{cvx}}$.

Approximate low-rank structure of the convex estimator $\bm{L}_{\mathsf{cvx}}$. Last but not least, Theorem 1 ensures that the convex estimate $\bm{L}_{\mathsf{cvx}}$ is nearly rank-$r$, so that a rank-$r$ approximation of $\bm{L}_{\mathsf{cvx}}$ is extremely accurate. In other words, the convex program automatically adapts to the true rank of $\bm{L}^{\star}$ without having any prior knowledge about $r$. As we shall see shortly, this is a crucial observation underlying the intimate connection between convex relaxation and a certain nonconvex approach.

Moving beyond the setting with $r,\kappa,\mu\asymp 1$, we have developed theoretical guarantees that allow $r,\kappa,\mu$ to grow with the problem dimension $n_{1},n_{2}$. The result is this.

Suppose that Assumptions 1-5 hold and that $n_{1}\geq n_{2}$. Take $\lambda=C_{\lambda}\sigma\sqrt{n_{1}p}$ and $\tau=C_{\tau}\sigma\sqrt{\log n_{2}}$ in (1.3) for some large enough constants $C_{\lambda},C_{\tau}>0$. Assume that

for some sufficiently large constant $C_{\mathsf{sample}}>0$ and some sufficiently small constants $c_{\mathsf{noise}},c_{\mathsf{outlier}}>0$. Then with probability exceeding $1-O(n_{2}^{-3})$, the following holds:

Any minimizer $(\bm{L}_{\mathsf{cvx}},\bm{S}_{\mathsf{cvx}})$ of the convex program (1.3) obeys

and the statistical guarantees (1.21) hold unchanged if $\bm{L}_{\mathsf{cvx}}$ is replaced by $\bm{L}_{\mathsf{cvx},r}$.

Similar to Theorem 1, our general theory (i.e. Theorem 2) provides the estimation error of the convex estimator $\bm{L}_{\mathsf{cvx}}$ in three different norms (i.e. the Euclidean, entrywise and operator norms), and reveals the near low-rankness of the convex estimator (cf. (1.22)) as well as the implicit regularization phenomenon (cf. (1.21b)).

Finally, we make note of several aspects of our general theory that call for further improvement. For instance, when there is no missing data and $n_{1}=n_{2}=n$, the rank $r$ of the unknown matrix $\bm{L}^{\star}$ needs to satisfy $r\lesssim\sqrt{n}$. On the positive side, our result allows $r$ to grow with the problem dimension $n$. However, prior results in the noiseless case [CLMW11, Li13] allow $r$ to grow almost linearly with $n$. This unsatisfactory aspect arises from the suboptimal analysis (in terms of the dependency on $r$) of a tightly related nonconvex estimation algorithm (to be elaborated on later), which, to the best of our knowledge, has not been resolved in the nonconvex low-rank matrix recovery literature [MWCC20, CLL20]. See Section 2 for more discussions about this point. Moreover, when $\bm{E}=\bm{0}$, it is known that $\rho_{s}$ can be as large as a constant even when the rank $r$ is allowed to grow with the dimension $n$ [Li13, CJSC13]. Our current theory, however, fails to cover the case with $\rho_{s}\asymp 1$ in the presence of noise. We demonstrate through numerical experiments that the dependence of $\rho_{\mathsf{s}}$ on $r$ might indeed by suboptimal in our current theory. More specifically, Figure 3 depicts the numerical Euclidean estimation errors w.r.t. the corruption probability $\rho_{\mathsf{s}}$ as we vary the rank while fixing the sampling ratio. It can be seen that the estimation error curves corresponding to different ranks align very well with each other, thus suggesting the capability of convex relaxation in tolerating a constant fraction $\rho_{\mathsf{s}}$ of outliers.

5 A peek at our technical approach

Before delving into the proof details, we immediately highlight our key technical ideas and novelties. For simplicity we assume $n_{1}=n_{2}=n$ throughout this section.

Instead of directly analyzing the convex program (1.3), we turn attention to a seemingly different, but in fact closely related, nonconvex program

It is worth emphasizing that our key idea — that is, bridging convex and nonconvex optimization — is drastically different from previous technical approaches for analyzing convex estimators (e.g. (1.3)). As it turns out, these prior approaches, which include constructing dual certificates and/or exploiting restricted strong convexity, have their own deficiencies in analyzing (1.3) and fall short of explaining the effectiveness of (1.3) in the random noise setting. For instance, constructing dual certificates in the noisy case is notoriously challenging given that we do not have closed-form expressions for the primal solutions (so that it is difficult to invoke the powerful dual construction strategies like the golfing scheme [Gro11] developed for the noiseless case). If we directly utilize the dual certificates constructed for the noiseless case, we would end up with an overly conservative bound like (1.4), which is exactly why the results in [ZLW+10, WL17] are sub-optimal. On the other hand, while it is viable to show certain strong convexity of (1.3) when restricted to some highly local sets and directions, it is unclear how (1.3) forces its solution to stay within the desired set and follow the desired directions, without adding further (and often unnecessary) constraints to (1.3).

It is worth noting that a similar connection between convex and nonconvex optimization has been pointed out by [CCF+20] towards understanding the power of convex relaxation for noisy matrix completion. Due to the absence of sparse outliers in the noisy matrix completion problem, the nonconvex loss function considered therein is smooth in nature, which greatly simplifies both the algorithmic and theoretical development. By contrast, the nonsmoothness inherent in (1.23) makes it particularly challenging to achieve the two desiderata mentioned above, namely, connecting the convex and nonconvex solutions and establishing the optimality of the nonconvex solution. In fact, to establish the connection between convex and nonconvex solutions, we put forward a novel two-step analysis strategy. Specifically, we first develop a crude upper bound on the Euclidean estimation error leveraging the idea of approximate dual certificates; see Theorem 3. While this crude upper bound is far from optimal, it serves as an important starting point towards formalizing the intimate relation between the convex solution $(\bm{L}_{\mathsf{cvx}},\bm{S}_{\mathsf{cvx}})$ and the nonconvex solution $(\bm{X},\bm{Y},\bm{S})$, since it is challenging to establish $\bm{X}\bm{Y}\approx\bm{L}_{\mathsf{cvx}}$ and $\bm{S}\approx\bm{S}_{\mathsf{cvx}}$ simultaneously without the aid of a crude bound. Second, in establishing the optimality of the nonconvex solution, the nonsmoothness nature of the nonconvex loss prevents us from applying the vanilla gradient descent scheme (as has been done in [CCF+20]). To address this issue, we develop an alternating minimization scheme — which alternates between gradient updates on $(\bm{X},\bm{Y})$ and minimization of $\bm{S}$ — aimed at minimizing the nonsmooth nonconvex loss function (1.23); see Algorithm 1 for details. As it turns out, such a simple algorithm allows us to track the proximity of the convex and nonconvex solutions and establish the optimality of the nonconvex solution all at once.

6 Random signs of outliers

The careful reader might wonder whether it is possible to remove the random sign assumption on $\bm{S}^{\star}$ (namely, Assumption 4) without compromising our statistical guarantees. After all, the results of [CSPW11, CLMW11, Li13] derived for the noise-free case do not rely on such a random sign assumption at all.Notably, in the noisy setting, prior theory [ZLW+10, WL17] also implicitly assumes this random sign condition, while [ANW12, KLT17] do not require this condition. Unfortunately, removal of such a condition might be problematic in general, as illustrated by the following example.

Suppose that (i) $n_{1}=n_{2}=n$, (ii) each non-zero entry of $\bm{S}^{\star}$ obeys $S_{ij}^{\star}=c_{0}\sigma$, (iii) $\rho_{\mathsf{s}}=c_{1}/\log n$ for some sufficiently small constant $c_{1}>0$, and (iv) there is no missing data (i.e. $p=1$). In such a scenario, the data matrix can be decomposed as

with high probability. Here the last step follows since $\widetilde{\bm{L}}^{\star}$ is of constant rank and condition number. This, however, leads to a lower bound on the estimation error

which can be $O(\sqrt{n}/\log n)$ times larger than the desired estimation error $O(\sigma\sqrt{n})$. Numerical experiments under the above setting (with $c_{0}=5$ and $c_{1}=1$) also suggest that (i) the estimation error under the fixed sign setting might be orderwise larger than that under the random sign setting; and (ii) under the fixed sign setting, the estimator (1.3) approximately recovers $\widetilde{\bm{L}}^{\star}$ instead of $\bm{L}^{\star}$; see Figure 5.

The take-away message is this: when the entries of $\bm{S}^{\star}$ are of non-random signs, it might sometimes be possible to decompose $\bm{S}^{\star}$ into (1) a low-rank bias component with a large Euclidean norm, and (2) a random fluctuation component whose typical size does not exceed that of $\bm{E}$. If this is the case, then the convex program (1.3) might mistakenly treat the bias component as a part of the low-rank matrix $\bm{L}^{\star}$, thus dramatically hampering its estimation accuracy.

Prior art

Principal component analysis (PCA) [Pea01, Jol11, FSZZ18] is one of the most widely used statistical methods for dimension reduction in data analysis. However, PCA is known to be quite sensitive to adversarial outliers — even a single corrupted data point can make PCA completely off. This motivated the investigation of robust PCA, which aims at making PCA robust to gross adversarial outliers. As formulated in [CLMW11, CSPW11], this is closely related to the problem of disentangling a low-rank matrix $\bm{L}^{\star}$ and a sparse outlier matrix $\bm{S}^{\star}$ (with unknown locations and magnitudes) from a superposition of them. Consequently, robust PCA can be viewed as an outlier-robust extension of the low-rank matrix estimation/completion tasks [CR09, KMO10, CLC19]. In a similar vein, robust PCA has also been extensively studied in the context of structured covariance estimation under approximate factor models [FFL08, FLM13, FWZ18, FWZ19], where the population covariance of certain random sample vectors is a mixture of a low-rank matrix and a sparse matrix, corresponding to the factor component and the idiosyncratic component, respectively.

Focusing on the convex relaxation approach, [CSPW11, CLMW11] started by considering the noiseless case with no missing data (i.e. $\bm{E}=\bm{0}$ and $p=1$) and demonstrated that, under mild conditions, convex relaxation succeeds in exactly decomposing both $\bm{L}^{\star}$ and $\bm{S}^{\star}$ from the data matrix $\bm{L}^{\star}+\bm{S}^{\star}$. More specifically, [CSPW11] adopted a deterministic model without assuming any probabilistic structure on the outlier matrix $\bm{S}^{\star}$. As shown in [CSPW11] and several subsequent work [CJSC13, HKZ11], convex relaxation is guaranteed to work as long as the fraction of outliers in each row/column does not exceed $O(1/r)$. In contrast, [CLMW11] proposed a random model by assuming that $\bm{S}^{\star}$ has random support (cf. Assumption 3); under this model, exact recovery is guaranteed even if a constant fraction of the entries of $\bm{S}^{\star}$ are nonzero with arbitrary magnitudes. Following the random location model proposed in [CLMW11], the paper [GWL+10] showed that, in the absence of noise, convex programming can provably tolerate a dominant fraction of outliers, provided that the signs of the nonzero entries of $\bm{S}^{\star}$ are randomly generated (cf. Assumption 4). Later, the papers [CJSC13, Li13] extended these results to the case when most entries of the matrix are unseen; even in the presence of highly incomplete data, convex relaxation still succeeds when a constant proportion of the observed entries are arbitrarily corrupted. It is worth noting that the results of [CJSC13] accommodated both models proposed in [CSPW11] and [CLMW11], while the results of [Li13] focused on the latter model.

The literature on robust PCA with not only sparse outliers but also dense noise — namely, when the measurements take the form $\bm{M}=\mathcal{P}_{\Omega_{\mathsf{obs}}}(\bm{L}^{\star}+\bm{S}^{\star}+\bm{E})$ — is relatively scarce. [ZLW+10, ANW12] were among the first to present a general theory for robust PCA with dense noise, which was further extended in [WL17, KLT17]. As we mentioned before, the first three [ZLW+10, ANW12, WL17] accommodated arbitrary noise with the last one [KLT17] focusing on the random noise. As we have discussed in Section 1.4, the statistical guarantees provided in these papers are highly suboptimal when it comes to the random noise setting considered herein. The paper [CC14] extended the robust PCA results to the case where the truth is not only low-rank but also of Hankel structure. The results therein, however, suffered from the same sub-optimality issue.

Moving beyond convex relaxation methods, another line of work proposed nonconvex approaches for robust PCA [NNS+14, GWL16, YPCC16, CGJ17, ZWG18, CCD+19, LMCC19, CCW19], largely motivated by the recent success of nonconvex methods in low-rank matrix factorization [CLC19, KMO10, CLS15, SL16, CC17, CW15, ZCL16, CC18, JNS13, NJS13, MWCC20, CCFM19, WGE17, WCCL16, CW18, ZL16, CDDD19]. Following the deterministic model of [CSPW11], the paper [NNS+14] proposed an alternating projection / minimization scheme to seek a low-rank and sparse decomposition of the observed data matrix. In the noiseless setting, i.e. $\bm{E}=\bm{0}$, this alternating minimization scheme provably disentangles the low-rank and sparse matrix from their superposition under mild conditions. In addition, [NNS+14] extended their result to the arbitrary noise case where the size of the noise is extremely small, namely, $\|\bm{E}\|_{\infty}\ll\sigma_{\min}/n$. When the noise $\{E_{ij}\}\sim\mathcal{N}(0,\sigma^{2})$, this is equivalent to the condition $\sigma\ll\sigma_{\min}/(n\sqrt{\log n})$. Comparing this with our noise condition $\sigma\ll\sigma_{\min}/(\sqrt{n\log n})$ (cf. (1.12)) when $r,\mu,\kappa\asymp 1$, one sees that our theoretical guarantees cover a wider range of noise levels. Similarly, [YPCC16] applied regularized gradient descent on a smooth nonconvex loss function which enjoys provable convergence guarantees to $(\bm{L}^{\star},\bm{S}^{\star})$ under the noiseless and partial observation setting. A recent paper [CCD+19] considered the nonsmooth nonconvex formulation for robust PCA and established rigorously the convergence of subgradient-type methods in the rank-1 setting, i.e. $r=1$. However, the extension to more general rank remains out of reach.

It is worth noting that noisy matrix completion problem [CP10, CCF+20] is subsumed as a special case by the model studied in this paper (namely, it is a special case with $\bm{S}^{\star}=\bm{0}$). Statistical optimality under the random noise setting (cf. Assumption 5) — including the convex relaxation approach [CCF+20, NW12, KLT11, Klo14] and the nonconvex approach [MWCC20, CLL20] — has been extensively studied. Focusing on arbitrary deterministic noise, [CP10] established the stability of the convex approach, whose resulting estimation error bound is similar to the one established for robust PCA with noise in [ZLW+10]) (see (1.4)). The paper [KS20] later confirmed that the estimation error bound established in [CP10] is the best one can hope for in the arbitrary noise setting for matrix completion, although it might be highly suboptimal if we restrict attention to random noise.

Finally, there is also a large literature considering robust PCA under different settings and/or from different perspectives. For instance, the computational efficiency in solving the convex optimization problem (1.3) and its variants has been studied in the optimization literature (e.g. [TY11, GMS13, SWZ14, MA18]). The problem has also been investigated under a streaming / online setting [GQV14, QV10, FXY13, ZLGV16, QVLH14, VN18]. These are beyond the scope of the current paper.

Architecture of the proof

In this section, we give an outline for proving Theorem 2. The proof of Theorem 1 follows immediately as it is a special case of Theorem 2. For simplicity of presentation, our proof sets $n_{1}=n_{2}=n$. It is straightforward to obtain the proof for the general rectangular case via minor modification.

The main ingredient of the proof lies in establishing an intimate link between convex and nonconvex optimization. Unless otherwise noted, we shall set the regularization parameters as

throughout. In addition, the soft thresholding operator at level $\tau$ is defined such that

For any matrix $\bm{X}$, the matrix $\mathcal{S}_{\tau}(\bm{X})$ is obtained by applying the soft thresholding operator $\mathcal{S}_{\tau}(\cdot)$ to each entry of $\bm{X}$ separately. Additionally, we define the true low-rank factors as follows

where $\bm{U}^{\star}\bm{\Sigma}^{\star}\bm{V}^{\star\top}$ is the SVD of the true low-rank matrix $\bm{L}^{\star}$.

We start by delivering a crude upper bound on the Euclidean estimation error, built upon the (approximate) duality certificate previously constructed in [CJSC13]. The proof is postponed to Appendix D.

Consider any given $\lambda>0$ and set $\tau\asymp\lambda\sqrt{({\log n})/{np}}$. Suppose that Assumptions 1-4 hold, and that

hold for some sufficiently large (resp. small) constant $C>0$ (resp. $c>0$). Then with probability at least $1-O(n^{-10})$, any minimizer $(\bm{L}_{\mathsf{cvx}},\bm{S}_{\mathsf{cvx}})$ of the convex program (1.3) satisfies

It is worth noting that the above theorem holds true for an arbitrary noise matrix $\bm{E}$. When specialized to the case with independent sub-Gaussian noise, this crude bound admits a simpler expression as follows.

Take $\lambda=C_{\lambda}\sigma\sqrt{np}$ and $\tau=C_{\tau}\sigma\sqrt{\log n}$ for some universal constant $C_{\lambda},C_{\tau}>0$. Under the assumptions of Theorem 3 and Assumption 5, we have — with probability exceeding $1-O(n^{-10})$ — that

This corollary follows immediately by combining Theorem 3 and Lemma 1 below. ∎

Suppose that Assumption 5 holds and that $n^{2}p>C_{1}n\log^{2}n$ for some sufficiently large constant $C_{1}>0$. Then with probability exceeding $1-O(n^{-10})$, one has

While the above results often lose a polynomial factor in $n$ vis-à-vis the optimal error bound, it serves as an important starting point that paves the way for subsequent analytical refinement.

2 Approximate stationary points of the nonconvex formulation

Instead of analyzing the convex estimator directly, we take a detour by considering the following nonconvex optimization problem

Here, $f\left(\bm{X},\bm{Y};\bm{S}\right)$ is a function of $\bm{X}$ and $\bm{Y}$ with $\bm{S}$ frozen, which contains the smooth component of the loss function $F(\bm{X},\bm{Y},\bm{S})$. As it turns out, the solution to convex relaxation (1.3) is exceedingly close to an estimate $(\bm{X},\bm{Y},\bm{S})$ obtained by a nonconvex algorithm aimed at solving (3.6) — to be detailed in Section 3.3. This fundamental connection between the two algorithmic paradigms provides a powerful framework that allows us to understand convex relaxation by studying nonconvex optimization.

In what follows, we set out to develop the afforementioned intimate connection. Before proceeding, we first state the following conditions concerned with the interplay between the noise size, the estimation accuracy of the nonconvex estimate $(\bm{X},\bm{Y},\bm{S})$, and the regularization parameters.

The regularization parameters $\lambda$ and $\tau\asymp\lambda\sqrt{({\log n})/{np}}$ satisfy

$\|\mathcal{P}_{\Omega_{\mathsf{obs}}}(\bm{E})\|<\lambda/16$ and $\|\mathcal{P}_{\Omega_{\mathsf{obs}}}(\bm{E})\|_{\infty}\leq\tau/4$;

$\|\bm{S}-\bm{S}^{\star}\|<\lambda/16$ and $\|\bm{X}\bm{Y}^{\top}-\bm{L}^{\star}\|_{\infty}\leq\tau/4$;

$\|\mathcal{P}_{\Omega_{\mathsf{obs}}}(\bm{X}\bm{Y}^{\top}-\bm{L}^{\star})-p(\bm{X}\bm{Y}^{\top}-\bm{L}^{\star})\|<\lambda/8$.

As an interpretation, the above condition says that: (1) the regularization parameters are not too small compared to the size of the noise, so as to ensure that we enforce a sufficiently large degree of regularization; (2) the estimate represented by the point $(\bm{X}\bm{Y}^{\top},\bm{S})$ is sufficiently close to the truth. At this point, whether this condition is meaningful or not remains far from clear; we shall return to justify its feasibility shortly.

Again, the validity of this condition will be discussed momentarily.

With the above conditions in place, we are ready to make precise the intimate link between convex relaxation and a candidate nonconvex solution. The proof is deferred to Appendix E.

Suppose that $n\geq\kappa$ and $\rho_{\mathsf{s}}\leq c/\kappa$ for some sufficiently small constant $c>0$. Assume that there exists a triple $(\bm{X},\bm{Y},\bm{S})$ such that

Further, assume that any singular value of $\bm{X}$ and $\bm{Y}$ lies in $[\sqrt{\sigma_{\min}/2},\sqrt{2\sigma_{\max}}]$. If the solution $(\bm{L}_{\mathsf{cvx}},\bm{S}_{\mathsf{cvx}})$ to the convex program (1.3) admits the following crude error bound

This theorem is a deterministic result, focusing on some sort of “approximate stationary points” of $F(\bm{X},\bm{Y},\bm{S})$. To interpret this, observe that in view of (3.7), one has $\nabla f\left(\bm{X},\bm{Y};\bm{S}\right)\approx\bm{0}$, and $\bm{S}$ minimizes $F(\bm{X},\bm{Y},\cdot)$ for any fixed $\bm{X}$ and $\bm{Y}$. If one can identify such an approximate stationary point that is sufficiently close to the truth (so that it satisfies Condition 1), then under mild conditions our theory asserts that

This would in turn formalize the intimate relation between the solution to convex relaxation and an approximate stationary point of the nonconvex formulation. The existence of such approximate stationary points will be verified shortly in Section 3.3.

The careful reader might immediately remark that this theorem does not say anything explicit about the minimizer of the nonconvex optimization problem (3.6); rather, it only pays attention to a special class of approximate stationary points of the nonconvex formulation. This arises mainly due to a technical consideration: it seems more difficult to analyze the nonconvex optimizer directly than to study certain approximate stationary points. Fortunately, our theorem indicates that any approximate stationary point obeying the above conditions serves as an extremely tight approximation of the convex estimate, and, therefore, it suffices to identify and analyze any such points.

3 Constructing an approximate stationary point via nonconvex algorithms

By virtue of Theorem 4, the key to understanding convex relaxation is to construct an approximate stationary point of the nonconvex problem (3.6) that enjoys desired statistical properties. For this purpose, we resort to the following iterative algorithm (Algorithm 1) to solve the nonconvex program (3.6).

In a nutshell, Algorithm 1 alternates between one iteration of gradient updates (w.r.t. the decision matrices $\bm{X}$ and $\bm{Y}$) and optimization of the non-smooth problem w.r.t. $\bm{S}$ (with $\bm{X}$ and $\bm{Y}$ frozen).Note that for any given $\bm{X}$ and $\bm{Y}$, the solution to $\text{minimize}_{\bm{S}}\ F(\bm{X},\bm{Y},\bm{S})$ is given precisely by $\mathcal{S}_{\tau}(\mathcal{P}_{\Omega_{\mathsf{obs}}}(\bm{M}-\bm{X}\bm{Y}^{\top}))$. For the sake of simplicity, we initialize this algorithm from the ground truth $(\bm{X}^{\star},\bm{Y}^{\star},\bm{S}^{\star})$, but our analysis framework might be extended to accommodate other more practical initialization (e.g. the one obtained by a spectral method [CCFM20]).

The following theorem makes precise the statistical guarantees of the above nonconvex optimization algorithm; the proof is deferred to Appendix F. Here and throughout, we define

where $\mathcal{O}^{r\times r}$ denotes the set of $r\times r$ orthonormal matrices.

Instate the assumptions of Theorem 2 and define

Take $t_{0}=n^{47}$ and $\eta\asymp 1/(n\kappa^{3}\sigma_{\max})$ in Algorithm 1. With probability at least $1-O(n^{-3})$, the iterates $\{(\bm{X}^{t},\bm{Y}^{t},\bm{S}^{t})\}_{0\leq t\leq t_{0}}$ of Algorithm 1 satisfy

In addition, with probability at least $1-O(n^{-3})$, one has

We shall also gather a few immediate consequences of Theorem 5 as follows, which contain basic properties that will be useful throughout.

Instate the assumptions of Theorem 5. Suppose that the sample size obeys $n^{2}p\gg\kappa^{4}\mu^{2}r^{2}n\log^{4}n$, the noise satisfies $\delta_{n}\ll 1/\sqrt{\kappa^{4}\mu r\log n}$, the outlier fraction satisfies $\rho_{\mathsf{s}}\ll 1/(\kappa^{3}\mu r\log n)$. With probability at least $1-O(n^{-3})$, the iterates of Algorithm 1 satisfy

4 Proof of Theorem 2

Theorem 5 and Corollary 2 have established appealing statistical performance of the nonconvex solution $(\bm{X}_{\mathsf{ncvx}},\bm{Y}_{\mathsf{ncvx}},\bm{S}_{\mathsf{ncvx}})$. To transfer this desired statistical property to that of $(\bm{L}_{\mathsf{cvx}},\bm{S}_{\mathsf{cvx}})$, it remains to show that the nonconvex estimator \big{(}\bm{X}_{\mathsf{ncvx}}\bm{Y}_{\mathsf{ncvx}}^{\top},\bm{S}_{\mathsf{ncvx}}\big{)} is extremely close to the convex estimator $(\bm{L}_{\mathsf{cvx}},\bm{S}_{\mathsf{cvx}})$. Towards this end, we intend to invoke Theorem 4; therefore, it boils down to verifying the conditions therein.

The small gradient condition (cf. (3.7)) holds automatically under (3.12).

By virtue of the spectral norm bound (3.11b), one has

as long as $\sigma\sqrt{\kappa n/p}\ll\sigma_{\min}$. This together with the Weyl inequality verifies the constraints on the singular values of $(\bm{X}_{\mathsf{ncvx}},\bm{Y}_{\mathsf{ncvx}})$.

The crude error bounds are valid in view of Theorem 3.

Regarding Condition 1 and Condition 2, Lemma 1 and standard inequalities about sub-Gaussian random variables imply that $\|\mathcal{P}_{\Omega_{\mathsf{obs}}}(\bm{E})\|<\lambda/16$ and $\|\mathcal{P}_{\Omega_{\mathsf{obs}}}(\bm{E})\|_{\infty}\leq\tau/4$. In addition, the bounds (3.11d) and (3.13b) ensure the second assumption $\|\bm{S}_{\mathsf{ncvx}}-\bm{S}^{\star}\|\leq\lambda/16$ and $\|\bm{X}\bm{Y}^{\top}-\bm{L}^{\star}\|_{\infty}\leq\tau/4$ in Condition 1. We are left with the last assumption in Condition 1 and Condition 2, which are guaranteed to hold in view of the following lemma (see Appendix C for the proof).

Instate the notations and assumptions of Theorem 2. Then with probability exceeding $1-O(n^{-10})$, we have

simultaneously for all $(\bm{X},\bm{Y})$ obeying

Here, $T$ denotes the tangent space of the set of rank-$r$ matrices at the point $\bm{X}\bm{Y}^{\top}$, and $C_{\infty}>0$ is an absolute constant.

Armed with the above conditions, we can readily invoke Theorem 4 to reach

with high probability. This taken collectively with Corollary 2 gives

Similar arguments lead to the advertised high-probability bounds

which establishes (1.22). In view of the triangle inequality, the properties (1.21) hold unchanged if $\bm{L}_{\mathsf{cvx}}$ is replaced by $\bm{L}_{\mathsf{cvx},r}$.

Discussion

This paper investigates the unreasonable effectiveness of convex programming in estimating an unknown low-rank matrix from grossly corrupted data. We develop an improved theory that confirms the optimality of convex relaxation in the presence of random noise, gross sparse outliers, and missing data. In particular, our results significantly improve upon the prior statistical guarantees [ZLW+10] under random noise, while further allowing for missing data. Our theoretical analysis is built upon an appealing connection between convex and nonconvex optimization, which has not been established previously.

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 grants W911NF-20-1-0097 and W911NF-18-1-0303, by the NSF grants CCF-1907661, IIS-1900140 and DMS-2014279, and by the Princeton SEAS innovation award. J. Fan is supported in part by the NSF grants DMS-1662139 and DMS-1712591, the ONR grant N00014-19-1-2120, and the NIH grant 2R01-GM072611-14.

Recall that $\Omega^{\star}$ is the support of the sparse component $\bm{S}^{*}$. In this section, we introduce an equivalent probabilistic model of $\Omega^{\star}$, which is more amenable to analysis and shall be assumed throughout the proof.

The original model. Recall from Assumption 3 the way we generate $\Omega^{\star}$ : (1) sample $\Omega_{\mathsf{obs}}$ from the i.i.d. Bernoulli model with parameter $p$; (2) for each $(i,j)\in\Omega_{\mathsf{obs}}$, let $(i,j)\in\Omega^{\star}$ independently with probability $\rho_{\mathsf{s}}$.

An equivalent model. The model involves over-sampling and rejection method: (1) sample $\Omega_{\mathsf{obs}}$ from the i.i.d. Bernoulli model with parameter $p$; (2) generate an augmented index set $\Omega_{\mathsf{aug}}\subseteq\Omega_{\mathsf{obs}}$ such that: for each $(i,j)\in\Omega_{\mathsf{obs}}$, we generate $(i,j)\in\Omega_{\mathsf{aug}}$ independently with probability $\rho_{\mathsf{aug}}$; (3) for any $(i,j)\in\Omega_{\mathsf{aug}}$, include $(i,j)$ in $\Omega^{\star}$ independently with probability $\rho_{\mathsf{s}}/\rho_{\mathsf{aug}}$.

It is straightforward to verify that the two models for $\Omega^{\star}$ are equivalent as long as $\rho_{\mathsf{s}}\leq\rho_{\mathsf{aug}}\leq 1$. Two important remarks are in order. First, by construction, we have $\Omega^{\star}\subseteq\Omega_{\mathsf{aug}}$. Second, the choice of $\rho_{\mathsf{aug}}$ can vary as needed as long as $\rho_{\mathsf{s}}\leq\rho_{\mathsf{aug}}\leq 1$.

Appendix B Preliminaries

This subsection collects several results that are useful throughout the proof. To begin with, the incoherence assumption (cf. Assumption 1) asserts that

where the first inequality comes from the elementary inequality $\|\bm{A}\bm{B}\|_{2,\infty}\leq\|\bm{A}\|_{2,\infty}\|\bm{B}\|$, and the last inequality is a consequence of the incoherence assumption as well as the fact that $\|(\bm{\Sigma}^{\star})^{1/2}\|=\|\bm{X}^{\star}\|$.

The next lemma is extensively used in the low-rank matrix completion literature.

Suppose that each $(i,j)$ is included in $\Omega_{0}\subseteq[n]\times[n]$ independently with probability $\rho_{0}$. Then with probability exceeding $1-O(n^{-10})$, one has

provided that $n^{2}\rho_{0}\gg\mu rn\log n$. Here, $T^{\star}$ denotes the tangent space of the set of rank-$r$ matrices at the point $\bm{L}^{\star}=\bm{X}^{\star}\bm{Y}^{\star\top}$.

In fact, the bound (B.2) uncovers certain near-isometry of the operator $\rho_{0}^{-1}\mathcal{P}_{\Omega_{0}}(\cdot)$ when restricted to the tangent space $T^{\star}$. This property is formalized in the following fact.

Suppose that $\|\mathcal{P}_{T^{\star}}-\rho_{0}^{-1}\mathcal{P}_{T^{\star}}\mathcal{P}_{\Omega_{0}}\mathcal{P}_{T^{\star}}\|\leq 1/2$. Then one has

The following corollary is an immediate consequence of Lemma 3 and Fact 1.

Suppose that $\rho_{\mathsf{s}}\leq\rho_{\mathsf{aug}}\leq 1/12$ and that $n^{2}p\rho_{\mathsf{aug}}\gg\mu rn\log n$. Then with probability at least $1-O(n^{-10})$, we have

Here, the second inequality arises from Lemma 3 and Fact 1 (by taking $\Omega_{0}=\Omega_{\mathsf{aug}}$ and $\rho_{0}=p\rho_{\mathsf{aug}}$). The proof is complete by recognizing the assumption $\rho_{\mathsf{aug}}\leq 1/12$. ∎

As it turns out, the near-isometry property of $\rho_{0}^{-1}\mathcal{P}_{\Omega_{0}}(\cdot)$ can be strengthened to a uniform version (uniform over a large collection of tangent spaces), as shown in the lemma below.

Suppose that each $(i,j)$ is included in $\Omega_{0}\subseteq[n]\times[n]$ independently with probability $\rho_{0}$, and that $n^{2}\rho_{0}\gg\mu rn\log n$. Then with probability at least $1-O(n^{-10})$,

holds simultaneously for all $(\bm{X},\bm{Y})$ obeying

Here, $c>0$ is some sufficiently small constant, and $T$ denotes the tangent space of the set of rank-$r$ matrices at the point $\bm{X}\bm{Y}^{\top}$.

Suppose that each $(i,j)$ is included in $\Omega_{0}\subseteq[n]\times[n]$ independently with probability $\rho_{0}$, and that $n^{2}\rho_{0}\gg n\log n$. Then there exists some absolute constant $C>0$ such that with probability at least $1-O(n^{-10})$,

holds simultaneously for all $\bm{A}$ and $\bm{B}$.

B.2 Proof of Lemma 4

The optimality condition of $(\bm{A},\bm{B})$ requires

see [CCF+20, Section C.3.1] for the justification of this identity. The proof then consists of two steps:

To see this, we can invoke the bound on $\alpha_{2}$ stated in [CCF+20, Appendix C.3.1] to yield

To this end, one starts with the following decomposition

In addition, the bound on $\alpha_{1}$ stated in [CCF+20, Appendix C.3.1] tells us that

Putting the above two bounds together, we conclude that

Appendix C Proof of Lemma 2

With Lemma 4 in place, we can immediately justify Lemma 2.

To begin with, the first two parts (3.16a) and (3.16b) are the same as [CCF+20, Lemma 4]. Hence, it suffices to verify the last one (3.16c). Recall from Appendix A that $\Omega^{\star}\subseteq\Omega_{\mathsf{aug}}$, where $\Omega_{\mathsf{aug}}$ is randomly sampled such that each $(i,j)$ is included in $\Omega_{\mathsf{aug}}$ independently with probability $p\rho_{\mathsf{aug}}$. Applying Lemma 4 on $\Omega_{\mathsf{aug}}$ finishes the proof, with the proviso that $\rho_{\mathsf{aug}}\asymp 1/\kappa^{2}$ and $\rho_{\mathsf{s}}\leq\rho_{\mathsf{aug}}$.

Appendix D Crude error bounds (Proof of Theorem 3)

Since $(\bm{L}_{\mathsf{cvx}},\bm{S}_{\mathsf{cvx}})$ is the minimizer of (1.3), it is self-evident that $\bm{S}_{\mathsf{cvx}}$ must be supported on $\Omega_{\mathsf{obs}}$. Then by construction, $\bm{\Lambda}_{\bm{S}},\bm{\Lambda}^{+}$ and $\bm{\Lambda}^{-}$ are all necessarily supported on $\Omega_{\mathsf{obs}}$, thus indicating that

Making use of this relation, we can continue the derivation (D.1) above to obtain

In the sequel, we shall control the three terms $\alpha_{1},\alpha_{2}$ and $\alpha_{3}$ separately.

Recognizing again that $\mathcal{P}_{\Omega_{\mathsf{obs}}}(\bm{L}^{\star}+\bm{S}^{\star}-\bm{M})=\mathcal{P}_{\Omega_{\mathsf{obs}}}(\bm{E})$, we can rearrange terms in (D.3) to derive

where we use the elementary inequality $a+b\leq\sqrt{2}\cdot\sqrt{a^{2}+b^{2}}$.

To relate $\alpha_{2}$ to $\alpha_{3}$, the following lemma plays a crucial role, whose proof is deferred to Appendix D.1.

Suppose that $\|\mathcal{P}_{\Omega^{\star}}\mathcal{P}_{T^{\star}}\|^{2}\leq p/8$ and that $\|\mathcal{P}_{T^{\star}}-p^{-1}\mathcal{P}_{T^{\star}}\mathcal{P}_{\Omega_{\mathsf{obs}}}\mathcal{P}_{T^{\star}}\|\leq 1/2$. Then for any pair $(\bm{A},\bm{B})$ of matrices, we have

The following lemma proves useful in linking $\alpha_{3}$ with $\alpha_{1}$, and we postpone the proof to Appendix D.2.

Here, the penultimate line results from the inequality (D.3) and last line follows from the same argument in obtaining (D.4).

This combined with (D.9) allows us to obtain

where we have used the elementary inequality $(a+b)^{2}\leq 2a^{2}+2b^{2}$. Recalling that $\tau=\lambda/\sqrt{np/\log n}$ and that $np\geq 1$, we arrive at

Taking the preceding bounds on $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$ collectively yields

Further, the elementary inequality $a^{2}+b^{2}\geq 2ab$ yields

We are left with proving that the conditions in Lemmas 6 and 7 hold with high probability. In view of Lemma 3 and Corollary 3, the conditions $\|\mathcal{P}_{T^{\star}}-p^{-1}\mathcal{P}_{T^{\star}}\mathcal{P}_{\Omega_{\mathsf{obs}}}\mathcal{P}_{T^{\star}}\|\leq 1/2$ and $\|\mathcal{P}_{\Omega^{\star}}\mathcal{P}_{T^{\star}}\|^{2}\leq p/8$ hold with high probability, provided that $n^{2}p\gg\mu rn\log n$ and $\rho_{\mathsf{s}}\leq 1/12$. In addition, Lemma 3 ensures that $\|\mathcal{P}_{T^{\star}}-p^{-1}(1-\rho_{\mathsf{s}})^{-1}\mathcal{P}_{T^{\star}}\mathcal{P}_{\Omega_{\mathsf{obs}}\setminus\Omega^{\star}}\mathcal{P}_{T^{\star}}\|\leq 1/2$ holds with high probability, with the proviso that $n^{2}p(1-\rho_{\mathsf{s}})\gg\mu rn\log n$, which holds true under the assumptions $\rho_{\mathsf{s}}\leq 1/12$ and $n^{2}p\gg\mu rn\log n$. Last but not least, the existence of the dual certificate $\bm{W}$ obeying (D.8) is guaranteed with high probability according to [CJSC13, Section III.D], under the conditions $\rho_{\mathsf{s}}\ll 1$ and $n^{2}p\gg\mu^{2}r^{2}n\log^{6}n$.Note that [CJSC13, Section III.D] requires $n^{2}p\gg\max\{\mu,\mu_{2}\}rn\log^{6}n$ under an additional incoherence condition $\|\bm{U}^{\star}\bm{V}^{\star\top}\|_{\infty}\leq\sqrt{\mu_{2}r/n^{2}}$. While we do not impose this extra condition, it is easily seen that $\|\bm{U}^{\star}\bm{V}^{\star\top}\|_{\infty}\leq\|\bm{U}^{\star}\|_{2,\infty}\|\bm{V}^{\star}\|_{2,\infty}\leq\mu r/n$ and hence $\mu_{2}\leq\mu^{2}r$.

D.1 Proof of Lemma 6

Here, the equality uses the fact $\Omega^{\star}\subseteq\Omega_{\mathsf{obs}}$, and the inequality holds because of the assumption $\|\mathcal{P}_{T^{\star}}-p^{-1}\mathcal{P}_{T^{\star}}\mathcal{P}_{\Omega_{\mathsf{obs}}}\mathcal{P}_{T^{\star}}\|\leq 1/2$ and Fact 1. Use $\Omega^{\star}\subseteq\Omega_{\mathsf{obs}}$ once again to obtain

Here, the last relation arises from the elementary inequality $ab\leq(a^{2}+b^{2})/2$ and the fact that $\left\|\mathcal{P}_{\Omega^{\star}}\mathcal{P}_{T^{\star}}\right\|\leq 1$. Combine the above two inequalities to obtain

as claimed, where we have used the assumption $\|\mathcal{P}_{\Omega^{\star}}\mathcal{P}_{T^{\star}}\|^{2}\leq p/8$ in the middle line and the fact $1/2\geq p/4$ in the last inequality.

D.2 Proof of Lemma 7

In view of the convexity of the nuclear norm $\|\cdot\|_{\ast}$, one has

where $\mathsf{sign}(\bm{S}^{\star})\coloneqq[\mathsf{sign}(S_{ij}^{\star})]_{1\leq i,j\leq n}$, and $\mathsf{sign}(\bm{S}^{\star})+\bm{G}_{2}$ is a sub-gradient of $\|\cdot\|_{1}$ at $\bm{S}^{\star}$. The first equality (i) holds by choosing $\bm{G}_{2}$ such that $-\langle\bm{G}_{2},\mathcal{P}_{\Omega_{\mathsf{obs}}}(\bm{H}_{\bm{L}})\rangle=\|\mathcal{P}_{\Omega_{\mathsf{obs}}\setminus\Omega^{\star}}(\bm{H}_{\bm{L}})\|_{1}$, and the last relation (ii) arises since $\mathsf{sign}(\bm{S}^{\star})$ is supported on $\Omega^{\star}\subseteq\Omega_{\mathsf{obs}}$. Combine the above two bounds to deduce that

In what follows, we shall lower bound the right-hand side of (D.11). To begin with, for $\theta_{1}$ we have

Here, the second identity uses the assumption (D.8c), the first inequality (i) uses the elementary inequality $|\langle\bm{A},\bm{B}\rangle|\leq\|\bm{A}\|_{\infty}\|\bm{B}\|_{1}$, and the last relation (ii) holds because of the assumption (D.8d). Substituting the above two bounds back into (D.11) gives

Take (D.13) and (D.14) collectively to yield

where the last relation is guaranteed by $np\gg 1$ and $\rho_{\mathsf{s}}\ll 1$. Recognizing that $\mathcal{P}_{\Omega_{\mathsf{obs}}\setminus\Omega^{\star}}(\bm{H}_{\bm{L}})=-\mathcal{P}_{\Omega_{\mathsf{obs}}\setminus\Omega^{\star}}(\bm{H}_{\bm{S}})$ finishes the proof.

Appendix E Equivalence between convex and nonconvex solutions (Proof of Theorem 4)

The goal of this section is to establish the intimate connection between the convex and nonconvex solutions (cf. Theorem 4). Before continuing, we remind the readers of the following notations:

$\bm{X}\bm{Y}^{\top}=\bm{U}\bm{\Sigma}\bm{V}^{\top}$: the rank-$r$ singular value decomposition of $\bm{X}\bm{Y}^{\top}$;

$T$: the tangent space of the set of rank-$r$ matrices at the estimate $\bm{X}\bm{Y}^{\top}$.

We begin with two useful lemmas which demonstrate that the point $(\bm{X}\bm{Y}^{\top},\bm{S})$ described in Theorem 4 satisfies approximate optimality conditions w.r.t. the convex program (1.3).

Instate the assumptions in Theorem 4. The triple $(\bm{X},\bm{Y},\bm{S})$ as stated in Theorem 4 satisfies

The proof can be straightforwardly adapted from [CCF+20, Claim 2] by replacing $\bm{E}$ therein with $\bm{E}+\bm{S}^{\star}-\bm{S}$. We omit it for the sake of brevity. ∎

The point $(\bm{X}\bm{Y}^{\top},\bm{S})$ as stated in Theorem 4 obeys

By definition, one has $\bm{S}=\mathcal{P}_{\Omega_{\mathsf{obs}}}[\mathcal{S}_{\tau}(\bm{M}-\bm{X}\bm{Y}^{\top})]$. Clearly, this is equivalent to saying that $\bm{S}$ is the unique minimizer of the following convex program

The claim of this lemma then follows from the optimality condition of this convex program (E.7). ∎

Additionally, in view of the crude error bound (3.8) and Condition 1, the matrix $\bm{\Delta}_{\bm{L}}$ (cf. (E.1)) obeys

E.2 Proof of Theorem 4

We now present the proof of Theorem 4, which consists of three main steps:

Showing that $(\bm{X}\bm{Y}^{\top},\bm{S})$ is not far from $(\bm{L}_{\mathsf{cvx}},\bm{S}_{\mathsf{cvx}})$ over $\Omega_{\mathsf{obs}}$, in the sense that $\mathcal{P}_{\Omega_{\mathsf{obs}}}(\bm{\Delta}_{\bm{L}}+\bm{\Delta}_{\bm{S}})\approx\bm{0}$;

Showing that $\bm{\Delta}_{\bm{L}}$ (resp. $\bm{\Delta}_{\bm{S}}$) is extremely small outside the tangent space $T$ (resp. the support $\Omega^{\star}$), and hence most of the energy of $\bm{\Delta}_{\bm{L}}$ (resp. $\bm{\Delta}_{\bm{S}}$) — if it is not vanishingly small — has to reside within $T$ (resp. $\Omega^{\star}$);

Showing that $\bm{\Delta}_{\bm{S}}\approx\bm{0}$ and $\bm{\Delta}_{\bm{L}}\approx\bm{0}$, with the assistance of the preceding two steps.

In what follows, we shall detail each of these steps.

Here, the equality arises from the relations $\bm{L}_{\mathsf{cvx}}=\bm{X}\bm{Y}^{\top}+\bm{\Delta}_{\bm{L}}$ and $\bm{S}_{\mathsf{cvx}}=\bm{S}+\bm{\Delta}_{\bm{S}}$. Expanding the squares and rearranging terms, we arrive at

Recall the definitions of $\bm{R}_{1}$ and $\bm{R}_{2}$ from Lemmas 8 and 9. We can then simplify the above inequality as

In the sequel, we develop bounds on $\theta_{1}$ and $\theta_{2}$.

With regards to $\theta_{1}$, one can further decompose it into

Similarly, one can decompose $\theta_{2}$ into

Combining (E.11), (E.12) and (E.13) yields

We begin by demonstrating that $\mathcal{P}_{T^{\perp}}(\bm{\Delta}_{\bm{L}})\approx\bm{0}$. From the inequality (E.14), we have

which demonstrates that the energy of $\bm{\Delta}_{\bm{L}}$ outside $T$ is extremely small.

where the relation holds since $\bm{\Delta}_{\bm{S}}$ is supported on $\Omega_{\mathsf{obs}}$. To facilitate the analysis of $\Omega_{\mathsf{obs}}\backslash\Omega^{\star}$, we introduce another index subset

The usefulness of $\Omega_{1}$ can be seen through the following claim, whose claim is postponed to the end of this section.

An immediate consequence of Claim 1 is that

which justifies our assertion that the energy of $\bm{\Delta}_{\bm{S}}$ outside $\Omega^{\star}$ is extremely small. Here, the last inequality arises from (E.15).

In view of (E.15) and the triangle inequality, we have

where the last step follows from (E.16) and (E.18). By Condition 2, we have

given that $\mathcal{P}_{T}(\bm{\Delta}_{\bm{L}})\in T$. The latter inequality combined with (E.15) and (E.16) further gives

Substituting the above bounds into (E.19) gives

provided that $n^{2}p\gg\kappa$. This combined with (E.16) allows one to control the size of $\bm{\Delta}_{\bm{L}}$:

In view of (E.15) and the fact that $\bm{\Delta}_{\bm{S}}$ is supported on $\Omega_{\mathsf{obs}}$, we have

E.2.4 Proof of Claim 1

where the second identity follows since $(\bm{L}_{\mathsf{cvx}},\bm{S}_{\mathsf{cvx}})=(\bm{X}\bm{Y}^{\top}+\bm{\Delta}_{\bm{L}},\bm{S}+\bm{\Delta}_{\bm{S}})$ is the optimizer of the convex program (1.3). These allow us to write

This characterization of $\bm{\Delta}_{\bm{S}}$ turns out to be crucial when establishing the inclusion $\Omega_{\mathsf{obs}}\backslash\Omega^{\star}\subseteq\Omega_{1}$. Towards this end, we need to introduce another index subset

As it turns out, the sets $\Omega,\Omega_{1}$ and $\Omega_{2}$ obey the following three conditions

Here, (i) follows since $\Omega\cup\Omega_{2}\subseteq\Omega^{\star}$, (ii) holds true since $\Omega_{2}\cap\Omega=\varnothing$, and (iii) results from the condition $\Omega_{\mathsf{obs}}\backslash\Omega\subseteq\Omega_{1}\cup\Omega_{2}$. It then boils down to proving each of the above three conditions.

The first one $\Omega_{2}\cap\Omega=\varnothing$ is straightforward to establish. Note that for any $(i,j)\in\Omega_{2}$, one must have \big{|}\big{(}\bm{M}-\bm{X}\bm{Y}^{\top}\big{)}_{ij}\big{|}\leq\tau and hence \big{[}\mathcal{S}_{\tau}(\bm{M}-\bm{X}\bm{Y}^{\top})\big{]}_{ij}=0, which means that $(i,j)\notin\Omega$. This proves the relation $\Omega_{2}\cap\Omega=\varnothing$.

Moving on to the second one $\Omega_{\mathsf{obs}}\backslash\Omega\subseteq\Omega_{1}\cup\Omega_{2}$, we prove this via contradiction. Suppose that this inclusion is false, i.e. there exits an index $(i,j)\in\Omega_{\mathsf{obs}}\backslash\Omega$ such that

Here, we have taken into account the fact that

To reach contradiction, we find it convenient to state the following simple fact.

Suppose that $|a|\leq\tau$ and that $\mathcal{S}_{\tau}(a+b)\neq 0$. Then

Given that $\mathcal{S}_{\tau}(a+b)\neq 0$, one necessarily has $|a+b|>\tau$. Without loss of generality, assume that $a+b>0$, which gives

This together with the fact $\tau\geq|a|$ yields $|\mathcal{S}_{\tau}(a+b)|=a+b-\tau\leq|b|+|a|-\tau$. ∎

With this fact in mind, we can deduce that

where (i) holds true since \big{[}\mathcal{S}_{\tau}(\bm{M}-\bm{X}\bm{Y}^{\top})\big{]}_{ij}=0 for any $(i,j)\in\Omega_{\mathsf{obs}}\backslash\Omega$, (ii) follows from Fact 2 (by taking $a=(\bm{M}-\bm{X}\bm{Y}^{\top})_{ij}$ and $b=\left[\bm{\Delta}_{\bm{S}}-\left(\bm{\Delta}_{\bm{L}}+\bm{\Delta}_{\bm{S}}\right)\right]_{ij}$), and (iii) is a consequence of (E.21) as well as the triangle inequality. The inequality (E.22), however, is clearly impossible. This establishes that $\Omega_{\mathsf{obs}}\backslash\Omega\subseteq\Omega_{1}\cup\Omega_{2}$.

We are left with the last one $\Omega\cup\Omega_{2}\subseteq\Omega^{\star}$, which is equivalent to saying $\Omega\subseteq\Omega^{\star}$ and $\Omega_{2}\subseteq\Omega^{\star}$. First, for any $(i,j)\in\Omega$, one has

Here, the last step comes from the triangle inequality and Condition 1. This reveals that $\Omega\subseteq\Omega^{\star}$. Similarly, for any $(i,j)\in\Omega_{2}$ we have

where we have used Condition 1, the bound (E.15), and the fact that $\tau\gg\sigma$. This demonstrates that $\Omega_{2}\subseteq\Omega^{\star}$. We have therefore justified that $\Omega\cup\Omega_{2}\subseteq\Omega^{\star}$.

Appendix F Analysis of the nonconvex procedure (Proof of Theorem 5)

This section is devoted to establishing Theorem 5. For notational convenience, we introduce

These allow us to express succinctly the rotation matrix $\bm{H}^{t}$ defined in (3.10) as

With the definitions of $\bm{F}^{t}$ and $\bm{H}^{t}$ in mind, it suffices to justify that: for all $0\leq t\leq t_{0}=n^{47}$, the following hypotheses

holds for all $1\leq t\leq t_{0}=n^{47}$.

Clearly, the bounds (3.11a), (3.11b), (3.11c), and (3.11d) in Theorem 5 follow immediately from (F.3a), (F.3b), (F.3c), and (F.3e), respectively. It remains to justify the small gradient bound (3.12) on the basis of (F.3) and (F.4), which is exactly the content of the following lemma.

Set $\lambda=C_{\lambda}\sigma\sqrt{np\log n}$ for some large constant $C_{\lambda}>0$. Suppose that $n^{2}p\gg\kappa^{3}\mu rn\log^{2}n$ and that the noise satisfies $\frac{\sigma}{\sigma_{\min}}\sqrt{\frac{n}{p}}\ll\frac{1}{\sqrt{\kappa^{4}\mu r\log n}}$. Take $\eta\asymp 1/(n\kappa^{3}\sigma_{\max})$. If the iterates satisfy (F.3) for all $0\leq t\leq t_{0}$ and (F.4) for all $1\leq t\leq t_{0}$, then with probability at least $1-O(n^{-50})$, one has

The remainder of this section is thus dedicated to showing that (F.3) and (F.4) hold for $\{(\bm{F}^{t},\bm{S}^{t})\}_{0\leq t\leq t_{0}}$, which we accomplish via mathematical induction. Throughout this section, we let $\bm{X}_{l,\cdot}$ denote the $l$th row of a matrix $\bm{X}$.

For each $1\leq l\leq n$, we define the following auxiliary loss functions

The above auxiliary loss function is obtained by dropping the randomness coming from the $l$th row of $\bm{M}$, which, as we shall see shortly, facilitates analysis in establishing the incoherence properties (F.3c). Similarly, we define for each $n+1\leq l\leq 2n$ that

Again, this auxiliary loss function is produced in a way that is independent from the $(l-n)$-th column of $\bm{M}$. In the above notation, $f^{\left(l\right)}\left(\bm{X},\bm{Y};\bm{S}\right)$ is a function of $\bm{X}$ and $\bm{Y}$ with $\bm{S}$ frozen.

For each $1\leq l\leq 2n$, we construct a sequence of leave-one-out iterates $\{\bm{F}^{t,(l)},\bm{S}^{t,(l)}\}_{t\geq 0}$ via Algorithm 2.

There are several features of the leave-one-out sequences that prove useful for our statistical analysis: (1) for the $l$th leave-one-out sequence, one can exploit the statistical independence to control the estimation error of $\bm{F}^{t,(l)}$ in the $l$th row; (2) the leave-one-out sequences and the original sequence $(\bm{F}^{t},\bm{S}^{t})$ are exceedingly close (since we have only discarded a small amount of information). These properties taken collectively allow us to control the estimation error of $\bm{F}^{t}$ in each row. To formalize these features, we make an additional set of induction hypotheses

Here, the rotation matrices $\bm{H}^{t,(l)}$ and $\bm{R}^{t,(l)}$ are defined respectively by

F.2 Key lemmas for establishing the induction hypotheses

This subsection establishes the induction hypotheses made in Appendix F.1, namely (F.3), (F.4) and (F.6). Before continuing, we find it convenient to introduce another function of $\bm{X}$ and $\bm{Y}$ (with $\bm{S}$ frozen) as follows

The difference between $f_{\mathsf{aug}}$ and $f$ lies in the following balancing term

that is, $f=f_{\mathsf{aug}}+f_{\mathsf{diff}}$.

The following four lemmas, which are inherited from [CCF+20] with little modification, are concerned with local strong convexity as well as the hypotheses (F.3a), (F.3b), (F.3d), (F.6b) and (F.3c).

Set $\lambda=C_{\lambda}\sigma\sqrt{np}$ for some large enough constant $C_{\lambda}>0$. Suppose that the sample size obeys $n^{2}p\gg\kappa\mu rn\log n$ and that the noise satisfies $\frac{\sigma}{\sigma_{\min}}\sqrt{\frac{n\log n}{p}}\ll 1$. Let the function $f_{\mathsf{aug}}$ be defined in (F.7). Then with probability at least $1-O(n^{-100})$,

Set $\lambda=C_{\lambda}\sigma\sqrt{np}$ for some large enough constant $C_{\lambda}>0$. Suppose that the sample size obeys $n^{2}p\gg\kappa\mu rn\log^{2}n$ and that the noise satisfies $\frac{\sigma}{\sigma_{\min}}\sqrt{\frac{n}{p}}\ll\frac{1}{\sqrt{\kappa^{4}\mu r\log n}}$. If the iterates satisfy (F.3) in the $t$th iteration, then with probability at least $1-O(n^{-100})$ ,

holds as long as $0<\eta\ll 1/(\kappa^{5/2}\sigma_{\max})$.

Set $\lambda=C_{\lambda}\sigma\sqrt{np}$ for some large enough constant $C_{\lambda}>0$. Suppose that the sample size obeys $n^{2}p\gg\kappa^{4}\mu^{2}r^{2}n\log^{2}n$ and that the noise satisfies $\frac{\sigma}{\sigma_{\min}}\sqrt{\frac{n}{p}}\ll\frac{1}{\sqrt{\kappa^{4}\log n}}$. If the iterates satisfy (F.3) in the $t$th iteration, then with probability at least $1-O(n^{-100})$, one has

Set $\lambda=C_{\lambda}\sigma\sqrt{np}$ for some large enough constant $C_{\lambda}>0$. Suppose that the sample size obeys $n^{2}p\gg\kappa^{2}\mu^{2}r^{2}n\log n$ and that the noise satisfies $\frac{\sigma}{\sigma_{\min}}\sqrt{\frac{n}{p}}\ll\frac{1}{\sqrt{\kappa^{2}\log n}}$. If the iterates satisfy (F.3) in the $t$th iteration, then with probability at least $1-O(n^{-100})$,

Set $\lambda=C_{\lambda}\sigma\sqrt{np}$ for some large enough constant $C_{\lambda}>0$. Suppose that the sample size obeys $n^{2}p\gg\kappa^{4}\mu^{2}r^{2}n\log^{3}n$ and that the noise satisfies $\frac{\sigma}{\sigma_{\min}}\sqrt{\frac{n}{p}}\ll\frac{1}{\sqrt{\kappa^{2}\log n}}$. If the iterates satisfy (F.3) in the $t$th iteration, then with probability at least $1-O(n^{-100})$,