High-dimensional principal component analysis with heterogeneous missingness
Ziwei Zhu, Tengyao Wang, Richard J. Samworth
TL;DR
This paper introduces primePCA, a new method for high-dimensional PCA with heterogeneous missing data, which outperforms existing estimators especially in realistic, noisy, and noiseless scenarios by iteratively imputing missing entries.
Contribution
The paper proposes primePCA, an iterative algorithm that effectively handles heterogeneous missingness in high-dimensional PCA, with theoretical guarantees and improved empirical performance.
Findings
primePCA converges to zero error in noiseless cases
Theoretical guarantees depend on average missingness properties
Numerical studies show strong performance across scenarios
Abstract
We study the problem of high-dimensional Principal Component Analysis (PCA) with missing observations. In simple, homogeneous missingness settings with a noise level of constant order, we show that an existing inverse-probability weighted (IPW) estimator of the leading principal components can (nearly) attain the minimax optimal rate of convergence. However, deeper investigation reveals both that, particularly in more realistic settings where the missingness mechanism is heterogeneous, the empirical performance of the IPW estimator can be unsatisfactory, and moreover that, in the noiseless case, it fails to provide exact recovery of the principal components. Our main contribution, then, is to introduce a new method for high-dimensional PCA, called `primePCA', that is designed to cope with situations where observations may be missing in a heterogeneous manner. Starting from the IPW…
| (H1) hardImpute | |||
|---|---|---|---|
| softImpute(oracle) | |||
| primePCA_init | |||
| primePCA | |||
| (H2) hardImpute | |||
| softImpute(oracle) | |||
| primePCA_init | |||
| primePCA | |||
| (H3) hardImpute | |||
| softImpute(oracle) | |||
| primePCA_init | |||
| primePCA | |||
| (H4) hardImpute | |||
| softImpute(oracle) | |||
| primePCA_init | |||
| primePCA |
| (H1) hardImpute | |||||
|---|---|---|---|---|---|
| softImpute(oracle) | |||||
| primePCA_init | |||||
| primePCA | |||||
| (H2) hardImpute | |||||
| softImpute(oracle) | |||||
| primePCA_init | |||||
| primePCA | |||||
| (H3) hardImpute | |||||
| softImpute(oracle) | |||||
| primePCA_init | |||||
| primePCA | |||||
| (H4) hardImpute | |||||
| softImpute(oracle) | |||||
| primePCA_init | |||||
| primePCA |
| ID | Title | Artist | Genre |
|---|---|---|---|
| 1 | Your Hand In Mine | Explosions In The Sky | Rock |
| 2 | All These Things That I’ve Done | The Killers | Rock |
| 3 | Lady Marmalade | Christina Aguilera / Lil’ Kim/ | Pop |
| Mya / Pink | |||
| 4 | Here It Goes Again | Ok Go | Rock |
| 5 | I Hate Pretending (Album Version) | Secret Machines | Rock |
| 6 | No Rain | Blind Melon | Rock |
| 7 | Comatose (Comes Alive Version) | Skillet | Rock |
| 8 | Life In Technicolor | Coldplay | Rock |
| 9 | New Soul | Yael Naïm | Pop |
| 10 | Blurry | Puddle Of Mudd | Rock |
| 11 | Give It Back | Polly Paulusma | Pop |
| 12 | Walking On The Moon | The Police | Rock |
| 13 | Face Down (Album Version) | The Red Jumpsuit Apparatus | Rock |
| 14 | Savior | Rise Against | Rock |
| 15 | Swing Swing | The All-American Rejects | Rock |
| 16 | Without Me | Eminem | Rap |
| 17 | Almaz | Randy Crawford | Pop |
| 18 | Hotel California | Eagles | Rock |
| 19 | Hey There Delilah | Plain White T’s | Rock |
| 20 | Revelry | Kings Of Leon | Rock |
| 21 | Undo | Björk | Rock |
| 22 | You’re The One | Dwight Yoakam | Country |
| Rock | Pop | Electronic | Rap | Country | RnB | Latin | Others | |
|---|---|---|---|---|---|---|---|---|
| Population | ||||||||
| Outliers |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
High-dimensional principal component analysis with heterogeneous missingness
Ziwei Zhu∗, Tengyao Wang*∗,†* and Richard J. Samworth*∗*
∗Statistical Laboratory, University of Cambridge
†Department of Statistical Science, University College London
Abstract
We study the problem of high-dimensional Principal Component Analysis (PCA) with missing observations. In simple, homogeneous missingness settings with a noise level of constant order, we show that an existing inverse-probability weighted (IPW) estimator of the leading principal components can (nearly) attain the minimax optimal rate of convergence. However, deeper investigation reveals both that, particularly in more realistic settings where the missingness mechanism is heterogeneous, the empirical performance of the IPW estimator can be unsatisfactory, and moreover that, in the noiseless case, it fails to provide exact recovery of the principal components. Our main contribution, then, is to introduce a new method for high-dimensional PCA, called primePCA, that is designed to cope with situations where observations may be missing in a heterogeneous manner. Starting from the IPW estimator, primePCA iteratively projects the observed entries of the data matrix onto the column space of our current estimate to impute the missing entries, and then updates our estimate by computing the leading right singular space of the imputed data matrix. It turns out that the interaction between the heterogeneity of missingness and the low-dimensional structure is crucial in determining the feasibility of the problem. We therefore introduce an incoherence condition on the principal components and prove that in the noiseless case, the error of primePCA converges to zero at a geometric rate when the signal strength is not too small. An important feature of our theoretical guarantees is that they depend on average, as opposed to worst-case, properties of the missingness mechanism. Our numerical studies on both simulated and real data reveal that primePCA exhibits very encouraging performance across a wide range of scenarios.
1 Introduction
One of the ironies of working with Big Data is that missing data play an ever more significant role, and often present serious difficulties for analysis. For instance, a common approach to handling missing data is to perform a so-called complete-case analysis (Little and Rubin, 2014), where we restrict attention to individuals in our study with no missing attributes. When relatively few features are recorded for each individual, one can frequently expect a sufficiently large proportion of complete cases that, under an appropriate missing at random hypothesis, a complete-case analysis may result in only a relatively small loss of efficiency. On the other hand, in high-dimensional regimes where there are many features of interest, there is often such a small proportion of complete cases that this approach becomes infeasible. As a very simple illustration of this phenomenon, imagine an data matrix in which each entry is missing independently with probability . When , a complete-case analysis would result in around 95 of the individuals (rows) being retained, but even when we reach , only around 5 of rows will have no missing entries.
The inadequacy of the complete-case approach in many applications has motivated numerous methodological developments in the field of missing data over the past 60 years or so, including imputation (Ford, 1983; Rubin, 2004), factored likelihood (Anderson, 1957) and Expectation-Maximisation approaches (Dempster, Laird and Rubin, 1977); see, e.g., Little and Rubin (2014) for an introduction to the area. Recent years have also witnessed increasing emphasis on understanding the performance of methods for dealing with missing data in a variety of high-dimensional problems, including sparse regression (Loh and Wainwright, 2012; Belloni, Rosenbaum and Tsybakov, 2017), classification (Cai and Zhang, 2018b), sparse principal component analysis (Elsener and van de Geer, 2018) and covariance and precision matrix estimation (Lounici, 2014; Loh and Tan, 2018).
In this paper, we study the effects of missing data in one of the canonical problems of high-dimensional data analysis, namely dimension reduction via Principal Component Analysis (PCA). This is closely related to the topic of matrix completion, which has received a great deal of attention in the literature over the last decade or so (e.g. Candès and Recht, 2009; Candès and Plan, 2010; Keshavan, Montanari and Oh, 2010; Mazumder, Hastie and Tibshirani, 2010; Koltchinskii, Lounici and Tsybakov, 2011; Candès et al., 2011; Negahban and Wainwright, 2012). There, the focus is typically on accurate recovery of the missing entries, subject to a low-rank assumption on the signal matrix; by contrast, our focus is on estimation of the principal eigenspaces. Previously proposed methods for low-dimensional PCA with missing data include non-linear iterative partial least squares (Wold and Lyttkens, 1969), iterative PCA (Kiers, 1997; Josse and Husson, 2012) and its regularised variant (Josse et al., 2009); see Dray and Josse (2015) for a nice survey and comparative study. More broadly, the R-miss-tastic website https://rmisstastic.netlify.com/ provides a valuable resource on methods for handling missing data.
The importance of the problem of high-dimensional PCA with missing data derives from its manifold applications. For instance, in many commercial settings, one may have a matrix of customers and products, with entries recording the number of purchases. Naturally, there will typically be a high proportion of missing entries. Nevertheless, PCA can be used to identify items that distinguish the preferences of customers particularly effectively, to make recommendations to users of products they might like and to summarise efficiently customers’ preferences. Later, we will illustrate such an application, on the Million Song Dataset, where we are able to identify particular songs that have substantial discriminatory power for users’ preferences as well as other interesting characteristics of the user database. Other potential application areas include health data, where one may seek features that best capture the variation in a population, and where the corresponding principal component scores may be used to cluster individuals into subgroups (that may, for instance, receive different treatment regimens).
To formalise the problem we consider, suppose that the (partially observed) matrix is of the form
[TABLE]
for independent random matrices and , where is a low-rank matrix and is a noise matrix with independent and identically distributed subgaussian entries having zero mean and unit variance. The low-rank property of is encoded through the assumption that it is generated via
[TABLE]
where has orthonormal columns and is a random matrix (with ) having independent and identically distributed rows with mean zero.
We are interested in estimating the column space of , denoted by , which is also the -dimensional leading eigenspace of . Cho, Kim and Rohe (2017) considered a different but related model where in (2) is deterministic, and is not necessarily centred, so that is the top right singular space of . (By contrast, in our setting, , so the mean structure is uninformative for recovering .) In the context of -homogeneous missingness, where each entry of is observed independently with probability , Cho, Kim and Rohe (2017) proposed to estimate by , where is a simple estimator formed as the top eigenvectors of an inverse-probability weighted (IPW) version of the sample covariance matrix (here, the weighting is designed to achieve approximate unbiasedness). Our first contribution, in Section 2, is to provide a detailed, finite-sample analysis of this estimator in the model given by (1) and (2), with a noise level of constant order. The differences between the settings necessitate completely different arguments, and reveal in particular a new phenomenon in the form of a phase transition in the attainable risk bound for the loss function, i.e. the Frobenius norm of the diagonal matrix of the sines of the principal angles between and . Moreover, we also provide a minimax lower bound in the case of estimating a single principal component, which reveals that this estimator achieves the minimax optimal rate up to a poly-logarithmic factor.
While this appears to be a very encouraging story for the IPW estimator, it turns out that it is really only the starting point for a more complete understanding of high-dimensional PCA. For instance, in the noiseless case, the IPW estimator fails to provide exact recovery of the principal components. Moreover, it is the norm rather than the exception in applications that the missingness mechanism is heterogeneous, in the sense that the probability of observing entries of varies (often significantly) across columns. For instance, in recommendation systems, some products will typically be more popular than others, and hence we observe more ratings in those columns. As another example, in meta-analyses of data from several studies, it is frequently the case that some covariates are common across all studies, while others appear only in a reduced proportion of them. In Section 2.2, we present an example to show that PCA algorithms can break down entirely for such heterogeneous missingness mechanisms when individual coordinates in may be large in absolute value. Intuitively, if we do not observe the interaction between the th and th columns of , then we cannot hope to estimate the th or th rows of , and this will cause substantial error if these rows of contain significant signal. This example illustrates that it is only possible to handle heterogeneous missingness in high-dimensional PCA with additional structure, and indicates that it is natural to assume incoherence among the entries of ; i.e. no single coordinate of is too large in absolute value.
Our main contribution, then, is to propose a new, iterative algorithm, called primePCA (short for projected refinement for imputation of missing entries in Principal Component Analysis), in Section 3, to estimate under this incoherence assumption, even with heterogeneous missingness. The initialiser for this algorithm is a modified version of the simple estimator discussed above, where the modification accounts for potential heterogeneity. Each iteration of primePCA projects the observed entries of onto the column space of the current estimate of to impute missing entries, and then updates our estimate of by computing the leading right singular space of the imputed data matrix. We show that in the noiseless setting, i.e., , primePCA achieves exact recovery of the principal eigenspaces (with a geometric convergence rate) when the initial estimator is close to the truth and a sufficiently large proportion of the data are observed. Moreover, we also provide a performance guarantee for the initial estimator, showing that it satisfies the desired requirement with high probability, conditional on any observed missingness pattern. Code for our algorithm is available in the R package primePCA (Zhu, Wang and Samworth, 2019).
To the best of our knowledge, primePCA is the first method for high-dimensional PCA that is designed to cope with settings where the missingness mechanism is heterogeneous. Indeed, the previously mentioned works on high-dimensional PCA and other high-dimensional statistical problems with missing data have either focused on a uniform missingness setting or have imposed a lower bound on entrywise observation probabilities, which reduces to this uniform case. A key contribution of our work is to account explicitly for the the effect of a heterogeneous missingness mechanism, where the estimation error depends on average entrywise missingness rather than worst-case missingness; see the discussions after Theorem 4 and Proposition 2 below. In Section 4, the empirical performance of primePCA is compared both with that of the initialiser, and with a popular method for matrix completion called softImpute (Mazumder, Hastie and Tibshirani, 2010; Hastie et al., 2015), which solves a nuclear-norm regularised optimisation problem, and which can be adapted to provide an estimate of . It turns out that in many settings, primePCA outperforms the softImpute algorithm, even when the latter is allowed access to the oracle choice of regularisation parameter for each dataset. Our analysis of the Million Song Dataset is given in Section 5. In Section 6, we illustrate how some of the ideas in this work may be applied to other high-dimensional statistical problems involving missing data. Proofs of our main results are deferred to Section 7; auxiliary results and their proofs are given in Section 8.
1.1 Notation
For a positive integer , we write . For and , we define \|\mathbf{v}\|_{p}:=\bigl{(}\sum_{j=1}^{d}|v_{j}|^{p}\bigr{)}^{1/p} and . We let denote the unit Euclidean sphere in . Given , we denote their Hamming distance by . Also, we write , , for the standard basis vector along the th coordinate axis and for the all-one vector in .
Given , we write for the diagonal matrix whose th diagonal entry is . We let denote the set of symmetric matrices in and let denote the set of matrices in with orthonormal columns. For a matrix , and , we write \|\mathbf{A}\|_{p}:=\bigl{(}\sum_{i,j}|A_{ij}|^{p}\bigr{)}^{1/p} if and for its entrywise norm, and for its -to- operator norm, where is viewed as a representation of a linear map from to . We provide special notation for the (Euclidean) operator norm and the Frobenius norm by writing and , respectively, and also write for the nuclear norm. If has the eigendecomposition for some and , we write for its th largest eigenvalue and abuse terminology slightly to refer to the leftmost columns of as the top eigenvectors of when . Also, we write . For any with singular value decomposition , where , and , we write for its th largest singular value and refer to the leftmost columns of (resp. ) as the top left (resp. right) singular vectors of . The Moore–Penrose pseudoinverse of is defined as . If , we write for the matrix obtained by extracting the rows of that are in .
For two matrices , we write if is positive semidefinite. The Kronecker product of two matrices and is defined as the block matrix
[TABLE]
When and , the Hadamard product of and , denoted , is defined such that for any and . Moreover, we say that is the Hadamard inverse of if .
If , then the matrix of principal angles between and is given by \Theta(\mathbf{U},\mathbf{V}):=\operatorname{diag}\bigl{(}\cos^{-1}(\sigma_{1}),\ldots,\allowbreak\cos^{-1}(\sigma_{K})\bigr{)}, where ; we let be defined entrywise. We define the loss function
[TABLE]
For a real-valued random variable and , we define its (Orlicz) -norm by
[TABLE]
For a random vector taking values in and , we define its (Orlicz) -norm by
[TABLE]
and define a version that is invariant to invertible affine transformations by
[TABLE]
We say that a -dimensional random vector is sub-Gaussian if . For two distributions and defined on the same measurable space and such that is absolutely continuous with respect to , the Kullback–Leibler divergence from to is given by
[TABLE]
Finally, for , we write if there exists a universal constant such that , and, where and may depend on an additional variable , say, we write if there exists , depending only on , such that .
2 The inverse-probability weighted estimator
For notational simplicity, we will write throughout the paper. Let denote the event that the th entry of is observed. We define the revelation matrix by and the partially observed data matrix
[TABLE]
In this section, we consider the simple case, where entries of the data matrix are observed independently and completely at random (i.e., independent of ) with -homogeneous missingness probability. Thus, for all , and and are independent for .
For , let and denote the th rows of and respectively, and define . Writing and for its Hadamard inverse, we have that under the -homogeneous missingness mechanism, \mathbf{P}=p^{2}\bigl{\{}\mathbf{1}_{d}\mathbf{1}_{d}^{\top}-(1-p^{-1})\mathbf{I}_{d}\bigr{\}} and \mathbf{W}=p^{-2}\bigl{\{}\mathbf{1}_{d}\mathbf{1}_{d}^{\top}-(1-p)\mathbf{I}_{d}\bigr{\}}. Following Cho, Kim and Rohe (2017), we consider the following weighted sample covariance matrix:
[TABLE]
The reason for including the weight is to ensure that , so that is an unbiased estimator of . Related ideas appear in the work of Cai and Zhang (2018a) on high-dimensional covariance matrix estimation with missing data. There, the authors propose a ‘generalised sample covariance matrix’, where the covariance between any two dimensions and is estimated using only the observations for which both dimensions and were observed. In practice, is typically unknown and needs to be estimated. It is thus natural to consider the following plug-in estimator :
[TABLE]
where \widehat{}\mathbf{W}=\widehat{p}^{-2}\bigl{\{}\mathbf{1}_{d}\mathbf{1}_{d}^{\top}-(1-\widehat{p})\mathbf{I}_{d}\bigr{\}} and denotes the proportion of observed entries in . We let denote the top eigenvectors of .
2.1 Theory for homogeneous missingness
In order to describe our theoretical performance guarantee for , we first list our conditions on the underlying data generating mechanism. We assume that is generated according to (1), (2) and (3), where:
- (A1)
, and are independent; 2. (A2)
has independent and identically distributed rows with and ; 3. (A3)
has independent and identically distributed entries with , and ; 4. (A4)
for all ; 5. (A5)
has independent entries.
In many places in this work, we will think intuitively of and as constants. In particular, if has multivariate normal rows and has normal entries, then we can simply take . For , under the same normality assumptions, we have , so this intuition amounts to thinking of the variance of each component of our data as being of constant order.
The theorem below gives bounds on the expected proximity between and .
Theorem 1**.**
Assume (A1)–(A5) and that , . Write . Then there exists a universal constant such that
[TABLE]
In particular, if , then there exists , depending only on and , such that
[TABLE]
When are regarded as constants, Theorem 2 below shows that (6) is the minimax rate up to logarithmic factors when . Note that the condition is reasonable given the scaling requirement for consistency of the empirical eigenvectors (Shen et al., 2016; Wang and Fan, 2017; Johnstone and Lu, 2009). Indeed, Theorem 5.1 in Shen et al. (2016) shows that when , the top eigenvector of the sample covariance matrix estimator is consistent if and only if . If we regard as the effective sample size in our missing data PCA problem, then it is a sensible analogy to assume here, which implies that the condition holds for large , up to poly-logarithmic factors.
As mentioned in the introduction, Cho, Kim and Rohe (2017) considered the different but related problem of singular space estimation in a model in which , where is a matrix of the form for a deterministic matrix , whose rows are not necessarily centred. In this setting is the leading -dimensional right singular space of , and the same estimator can be applied. An important distinction is that, when the rows of are not centred and the entries of are of comparable magnitude, is of order , so when is regarded as a constant, it is natural to think of the singular values of as also being of order . Indeed, this is assumed in Cho, Kim and Rohe (2017). On the other hand, in our model, where the rows of have mean zero, assuming that the eigenvalues are of order would amount to an extremely strong requirement, essentially restricting attention to very highly spiked covariance matrices. Removing this condition requires completely different arguments. Moreover, (6) reveals an interesting phase transition phenomenon that has not been observed previously in the literature. Specifically, if the signal strength is large enough that , then we should regard as the effective sample size, as might intuitively be expected. On the other hand, if , then the estimation problem is considerably more difficult and the effective sample size is of order . In fact, by inspecting the proof of Theorem 1, we see that in the high signal case, it is the difficulty of estimating the diagonal entries of that drives the rate, while when the signal strength is low, the bottleneck is the challenge of estimating the off-diagonal entries. By comparing (6) with the minimax lower bound result in Theorem 2 below, we see that this phase transition phenomenon is an inherent feature of this estimation problem rather than an artefact of the proof techniques we used to derived the upper bound.
In order to state our minimax lower bound, we let denote the class of distributions of pairs satisfying (A1), (A2), (A3) and (A5) with . Since we are now working with vectors instead of matrices, we write in place of .
Theorem 2**.**
There exists a universal constant such that
[TABLE]
where the infimum is taken over all estimators of .
Theorem 2 reveals that in Theorem 1 achieves the minimax optimal rate of estimation up to a logarithmic factor when and are regarded as constants and .
2.2 General observation mechanism
Although Section 2.1 may appear to indicate that the problem of high-dimensional PCA with missing entries is essentially solved, the aim of this subsection is to show that the situation changes dramatically once the data can be missing heterogeneously.
To this end, consider the following example. Suppose that is equal to or with equal probability, so that
[TABLE]
In other words, for each , we observe precisely one of the first two entries of , together with all of the remaining entries. Let , where , and , where . Suppose that and let , and similarly assume that and set . Then and are identically distributed. However, the leading eigenvectors of and are respectively and , which are orthogonal! Thus, it is impossible to simultaneously estimate consistently the leading eigenvector of both and from our observations. We note that it is the disproportionate weight of the first two coordinates in the leading eigenvector, combined with the failure to observe simultaneously the first two entries in the data, that makes the estimation problem intractable in this example.
The understanding derived from this example motivates us to assume that satisfies the incoherence condition for some . The intuition here is that the maximally incoherent case is where each column of is a unit vector proportional to a vector whose entries are either or , in which case . Our condition asks for the columns of to have an incoherence of the same order as this maximally incoherent case. Similar conditions have been invoked in the literature on matrix completion (e.g., Candès and Plan, 2010; Keshavan, Montanari and Oh, 2010), but for a different reason. There, the purpose is to ensure that the true right singular space is not too closely aligned with the standard basis, which allows the missing entries of the matrix to be inferred from relatively few observations. In our case, the incoherence condition ensures that significant estimation error in a few components of the leading eigenvectors does not affect the overall statistical performance too much. Thus, it rules out examples such as the one described above, where heavy corruption in only a few entries spoils any chance of consistent estimation.
3 Our new algorithm, primePCA
We are now in a position to introduce and analyse our iterative algorithm primePCA to estimate principal eigenspaces of the covariance matrix . The basic idea is to iteratively refine a current (input) iterate by first imputing the missing entries of the data matrix using the current estimate of , and then applying a singular value decomposition (SVD) to the completed data matrix. More precisely, for , we let denote the indices for which the corresponding entry of is observed, and regress the observed data on to obtain an estimate of the th row of . This is natural in view of the data generating mechanism . We then use to impute the missing values , retain the original observed entries as , and set our next (output) iterate to be the top right singular vectors of the imputed matrix . To motivate this final choice, observe that when , we have ; we therefore have the SVD , where and is diagonal with positive diagonal entries. This means that , so the column spaces of and coincide. For convenience, pseudocode of a single iteration of refinement in this algorithm is given in Algorithm 1.
We now seek to provide formal justification for Algorithm 1. For any , we let be an SVD of and let . The two-to-infinity distance between and is then defined to be
[TABLE]
Note that the definition of does not depend on our choice of SVD and that for any . The following proposition considers the noiseless setting , and provides conditions under which, for any estimator that is close to , a single iteration of refinement in Algorithm 1 contracts the two-to-infinity distance between their column spaces. In a slight abuse of notation, we write .
Proposition 1**.**
Let \widehat{}\mathbf{V}_{K}^{(\mathrm{out})}:=\emph{{refine}}\bigl{(}K,\widehat{}\mathbf{V}_{K}^{(\mathrm{in})},{\boldsymbol{\Omega}},{\mathbf{Y}}_{{\boldsymbol{\Omega}}}\bigr{)} as in Algorithm 1 and let . We assume that , that , and write the SVD of as . Suppose that , and that both and hold for some . Then there exist , depending only on and , such that whenever
- (i)
, 2. (ii)
,
we have that
[TABLE]
In practice, in cases where either of the two conditions on or is not satisfied, we first perform a screening step that restricts attention to a set of row indices for which the data contain sufficient information to estimate the principal components. This screening step is explicitly accounted for in Algorithm 2 below, as well as in the theory that justifies it.
Algorithm 2 provides pseudocode for the iterative primePCA algorithm, given an initial estimator . The iterations continue until either we hit the convergence threshold or the maximum iteration number . Theorem 3 below guarantees that, in the noiseless setting of Proposition 1, the primePCA estimator converges to at a geometric rate.
Theorem 3**.**
For , let be the iterate of Algorithm 2 with input , , , , , and . Write and let
[TABLE]
where . Let be an SVD of . Suppose that both and hold, for some . Let
[TABLE]
and assume that . Then there exist , depending only on , and , such that whenever
- (i)
, 2. (ii)
,
we have that for every ,
[TABLE]
3.1 Initialisation
Theorem 3 provides a general guarantee on the performance of primePCA, but relies on finding an initial estimator that is sufficiently close to the truth . The aim of this subsection, then, is to propose a simple initialiser and show that it satisfies the requirement of Theorem 3 with high probability, conditional on the missingness pattern.
Consider the following modified weighted sample covariance matrix
[TABLE]
where for any ,
[TABLE]
Here, the matrix replaces in (4) because we no longer wish to assume homogeneous missingness. We take as our initial estimator of the matrix of top eigenvectors of , denoted . Theorem 4 below studies the performance of this initialiser, in terms of its two-to-infinity norm error, as required for application in Theorem 3. We write and for probabilities and expectations conditional on .
Theorem 4**.**
Assume (A1)–(A4) and that . Suppose further that , that for all and let . Then there exist , depending only on and , such that for every , if
[TABLE]
then
[TABLE]
As a consequence, writing
[TABLE]
where is as in Theorem 3, we have that
[TABLE]
The first part of Theorem 4 provides a general probabilistic bound for , after conditioning on the missingness pattern. This allows us, in the second part, to provide a guarantee on the probability with which is a good enough initialiser for Theorem 3 to apply. For intuition regarding , consider the -homogenous missingness setting. In that case, by Lemma 6, typical realisations of have and when , so in the spiked model where and are both of order , we expect to be small.
One of the attractions of our analysis is the fact that we are able to provide bounds that only depend on entrywise missingness probabilities in an average sense, as opposed to worst-case missingness probabilities. The refinements conferred by such bounds are particularly important when the missingness mechanism is heterogeneous, as typically encountered in practice. The averaging of missingness probabilities can be partially seen in Theorem 4, since and depend only on the and norms of each row of , but is even more evident in the proposition below, which gives a probabilistic bound on the original distance between and .
Proposition 2**.**
Assume the same conditions as in Theorem 4. Then there exists a universal constant such that for any , if
[TABLE]
*then *
[TABLE]
In this bound, then, we see that only depends on through the entrywise and norms of the whole matrix. Lemma 6 provides probabilistic control of these norms under the -homogeneous missingness mechanism. In general, if the rows of are independent and identically distributed, but different covariates are missing with different probabilities, then entries of will concentrate around the reciprocals of the simultaneous observation probabilities of pairs of covariates. As such, for a typical realisation of , our bound in Proposition 2 depends only on the harmonic averages of these simultaneous observation probabilities and their squares. Such an averaging effect ensures that our method is effective in a much wider range of heterogeneous settings than previously allowed in the literature.
4 Simulation studies
In this section, we assess the empirical performance of primePCA as proposed in Algorithm 2, with initialiser from Section 3.1, and denote the output of this algorithm by . We generate observations according to the model described in (1), (2) and (3) where the rows of the matrix are independent random vectors, for some . We further generate the observation indicator matrix , independently of and , and investigate the following four missingness mechanisms that represent different levels of heterogeneity:
- (H1)
Homogeneous: for all ; 2. (H2)
Mildly heterogeneous: for , where and independently; 3. (H3)
Highly heterogeneous columns: for and all odd and for and all even . 4. (H4)
Highly heterogeneous rows: for and all odd and for and all even .
In Sections 4.1, 4.2 and 4.3 below, we investigate primePCA in noiseless, noisy and misspecified settings respectively. In all cases, the average statistical error was estimated from 100 Monte Carlo repetitions of the experiment. For comparison, we also studied the softImpute algorithm (Mazumder, Hastie and Tibshirani, 2010; Hastie et al., 2015), which is considered to be state-of-the-art for matrix completion (Chi, Lu and Chen, 2018). This algorithm imputes the missing entries of by solving the following nuclear-norm-regularised optimisation problem:
[TABLE]
where is to be chosen by the practitioner. The softImpute estimator of is then given by the top right singular space of .
Figure 1 presents Monte Carlo estimates of for different choices of in two different settings. The first uses the noiseless set-up of Section 4.1, together with missingness mechanism (H1); the second uses the noisy setting of Section 4.2 with parameter and missingness mechanism (H2). We see that the error barely changes when varies within ; very similar plots were obtained for different data generation and missingness mechanisms, though we omit these for brevity. For definiteness, we therefore fixed throughout our simulation study.
4.1 Noiseless case
In the noiseless setting, we set , and also fix , , and . We set
[TABLE]
In Figure 2, we present the logarithm of the estimated average loss of primePCA and softImpute under (H1), (H2), (H3) and (H4). We set the range of -axis to be the same for each method to facilitate straightforward comparison. We see that the statistical error of primePCA decreases geometrically as the number of iterations increases, which confirms the conclusion of Theorem 3 in this noiseless setting. Moreover, after a moderate number of iterations, its performance is a substantial improvement on that of the softImpute algorithm, even if this latter algorithm is given access to an oracle choice of the regularisation parameter . The high statistical error of softImpute in these settings can be partly explained by the default value of the tuning parameter thresh in the softImpute package in R, namely , which corresponds to the red curve in the right-hand panels of Figure 2. By reducing the values of thresh to and , corresponding to the green and blue curves in Figure 2 respectively, we were able to improve the performance of softImpute to some extent, though the statistical error is sensitive to the choice of the regularisation parameter . Moreover, even with the optimal choice of , it is not competitive with primePCA (which is also considerably faster to compute, even with iterations).
4.2 Noisy case
Here, we generate the rows of as independent random vectors, independent of all other data. We maintain the same choices of , , and as in Section 4.1, set and vary to achieve different signal-to-noise ratios. In particular, defining , the choices correspond to the low, medium and high signal-to-noise ratios , respectively. For an additional comparision, we also considered a variant of the softImpute algorithm called hardImpute (Mazumder, Hastie and Tibshirani, 2010), which retains only a fixed number of top singular values in each iteration of matrix imputation; this can be achieved by setting the argument in the softImpute function to be [math].
We remark that in general, the choice of for softImpute is more challenging than in many regularised -estimation contexts, because in our setting we have no response variable, so cross-validation techniques are less readily available. For our comparisons, therefore, we gave the softImpute algorithm a particularly strong form of oracle choice of , namely where was chosen for each individual repetition of the experiment, so as to minimise the loss function. Naturally, such a choice is not available to the practitioner. Moreover, in order to ensure the range of was wide enough to include the best softImpute solution, we set the argument rank.max in that algorithm to be .
In Table 1, we report the statistical error of primePCA after iterations of refinement, together with the corresponding statistical errors of our initial estimator primePCA_init and those of softImpute(oracle) and hardImpute. Remarkably, primePCA exhibits stronger performance than these other methods across each of the signal-to-noise ratio regimes and different missingness mechanisms. We also remark that hardImpute is inaccurate and unstable, because it might converge to the local optimum that is far from the truth.
4.3 Near low-rank case
Here, we set , , , , and fixed once for all experiments to be the top eigenvectors of one realisation111In R, we set the random seed to be before generating . of the sample covariance matrix of independent random vectors. Here , and we again generated the rows of as independent random vectors. Table 2 reports the average loss of estimating the top eigenvectors of , where varies from to . Interestingly, even in this misspecified setting, primePCA is competitive with the oracle version of softImpute.
5 Real data analysis: Million Song Dataset
We apply primePCA to a subset of the Million Song Dataset222https://www.kaggle.com/c/msdchallenge/data to analyse music preferences. The original data can be expressed as a matrix with users (rows) and songs (columns), with entries representing the number of times a song was played by a particular user. The proportion of non-missing entries in the matrix is . Since the matrix is very sparse, and since most songs have very few listeners, we enhance the signal-to-noise ratio by restricting our attention to songs that have at least 100 listeners ( songs in total). This improves the proportion of non-missing entries to . Further summary information about the filtered data is provided below:
Quantiles of non-missing matrix entry values:
[TABLE] 2. 2.
Quantiles of the number of listeners for each song:
[TABLE] 3. 3.
Quantiles of the total play counts of each user:
[TABLE]
Moreover, from the first numbered point above, we see that the distribution of play counts has an extremely heavy tail. To guard against excessive influence from the outliers, we discretise the play counts into five interest levels as follows:
[TABLE]
We are now in a position to analyse the data using primePCA. For and , let denote the level of interest of user in song , let and let . Our initial goal is to assess the top eigenvalues of to see if there is low-rank signal in . To this end, we first apply Algorithm 2 to obtain ; next, for each , we run Steps 2–5 of Algorithm 1 to obtain the estimated principal score , so that we can approximate by . This allows us to estimate by . Figure 3 displays the top eigenvalues of , which exhibit a fairly rapid decay, thereby providing evidence for the existence of low-rank signal in .
In the left panel of Figure 4, we present the estimate of the top two eigenvectors of the covariance matrix , with colours indicating the genre of the song. The outliers in the -axis of this plot are particularly interesting: they reveal songs that polarise opinion among users (see Table 3) and that best capture variation in individuals’ preferences for types of music measured by the first principal component. It is notable that Rock songs are overrepresented among the outliers (see Table 4), relative to, say, Country songs. Users who express a preference for particular songs are also more likely to enjoy songs that are nearby in the plot. Such information is therefore potentially commercially valuable, both as an efficient means of gauging users’ preferences, and for providing recommendations.
The right panel of Figure 4 presents the principal scores of the users, with frequent users (whose total song plays are in the top 10% of all users) in red and occasional users in blue. This plot reveals, for instance, that the second principal component is well aligned with general interest in the website. Returning to the left plot, we can now interpret a positive -coordinate for a particular song (which is the case for the large majority of songs) as being associated with an overall interest in the music provided by the site.
6 Discussion
Heterogeneous missingness is ubiquitous in contemporary, large-scale data sets, yet we currently understand very little about how existing procedures perform or should be adapted to cope with the challenges this presents. Here we attempt to extract the lessons learned from this study of high-dimensional PCA, in order to see how related ideas may be relevant in other statistical problems where one wishes to recover low-dimensional structure with data corrupted in a heterogeneous manner.
A key insight, as gleaned from Section 2.2, is that the way in which the heterogeneity interacts with the underlying structure of interest is crucial. In the worst case, the missingness may be constructed to conceal precisely the structure one seeks to uncover, thereby rendering the problem infeasible by any method. The only hope, then, in terms of providing theoretical guarantees, is to rule out such an adversarial interaction. This was achieved via our incoherence condition in Section 3, and we look forward to seeing how the relevant interactions between structure and heterogeneity can be controlled in other statistical problems such as those mentioned in the introduction. For instance, in sparse linear regression, one would anticipate that missingness of covariates with strong signal would be much more harmful than corresponding missingness for noise variables.
Our study also contributes to the broader understanding of the uses and limitations of spectral methods for estimating hidden low-dimensional structures in high-dimensional problems. We have seen that the IPW estimator is both methodologically simple and achieves near-minimax optimality when the noise level is of constant order. Similar results have been obtained for spectral clustering for network community detection in stochastic block models (Rohe et al., 2011) and in low-rank-plus-sparse matrix estimation problems (Fan, Liao and Mincheva, 2013). On the other hand, the IPW estimator fails to provide exact recovery of the principal components in the noiseless setting. In these other aforementioned problems, it has also been observed that refinement of an initial spectral estimator can enhance performance, particularly in high signal-to-noise ratio regimes (Gao et al., 2016; Zhang, Cai and Wu, 2018), as we were able to show for our primePCA algorithm. This suggests that such a refinement has the potential to confer a sharper dependence of the statistical error rate on the signal-to-noise ratio compared with a vanilla spectral algorithm, and understanding this phenomenon in greater detail provides another interesting avenue for future research.
7 Proofs of main results
We define two linear maps , such that for any , we have and . In other words, and correspond to the diagonal and off-diagonal parts of respectively.
Proof of Theorem 1.
Since , we have that
[TABLE]
Moreover, since by Lemma 1, it follows from van der Vaart and Wellner (1996, Lemma 2.2.2) that there exist a universal constant such that333In van der Vaart and Wellner (1996), the -norm of a random variable is defined slightly differently as . It can be shown (Vershynin, 2012, Lemma 5.5) that these two norms are equivalent.
[TABLE]
Recall that denotes the th row of . Define and . We have the following decomposition of :
[TABLE]
We regard as a perturbed version of . Applying Yu, Wang and Samworth (2015, Theorem 2), we have
[TABLE]
We will control the expectation of the three terms on the right-hand side of (7) separately. Define . For notational simplicity, we write and respectively for the probability and expectation conditional on . Also, let and (if , then ). For the first term, we apply a symmetrisation argument. Let denote copies of that are independent of , let be independent Rademacher random variables that are independent of and write for expectation conditional on . Then by Jensen’s inequality,
[TABLE]
Since , we have that
[TABLE]
Writing , we then have
[TABLE]
Notice that
[TABLE]
Therefore,
[TABLE]
Now, observe that , so for ,
[TABLE]
By the Cauchy–Schwarz inequality, we therefore have that
[TABLE]
where is a universal constant, the second inequality uses (9) and (10) and the penultimate bound uses Stirling’s inequality.
Let and \sigma^{2}:=C^{\prime}Mn^{-1}d\log d\sum_{i=1}^{n}\bigl{\{}R\tau^{2}\widehat{p}_{i}^{(3)}+\widehat{p}_{i}^{(2)}M\log d\bigr{\}}. Then by Tropp (2012, Theorem 6.2), we obtain that
[TABLE]
Consequently, for , we have
[TABLE]
Given (7), integrating the left-hand side of the above inequality over yields
[TABLE]
where the first inequality uses Jensen’s inequality and the second inequality uses Lemma 4.
For the second sum on the right-hand side of (7), we have by van der Vaart and Wellner (1996, Lemma 2.2.2) again that
[TABLE]
where the final inequality uses Lemma 2 and the fact that . Now by the Cauchy–Schwarz inequality,
[TABLE]
which is dominated by the bound in (7).
Finally, for the third term on the right-hand side of (7), we have by the Cauchy–Schwarz inequality again that
[TABLE]
which is also dominated by the bound in (7). Substituting (7), (7) and (15) into (7) establishes (5). If we regard and as constants and if , then the second term in the bracket of the right-hand side of (5) is dominated by the first term, and claim (6) follows immediately. ∎
Proof of Theorem 2.
Without loss of generality, we may assume that and that is even, and write for some . By the Gilbert–Varshamov lemma (see, e.g. Massart, 2007, Lemma 4.7), there exist such that and for any distinct pair of vectors . Let be a real number to be specified later. To each , we can associate a distribution such that is a random vector ( random matrix) with independent entries, is an random matrix with independent entries, and
[TABLE]
Fixing distinct , we write and and let and denote respectively the marginal distribution of under and . Define and also set and . Then
[TABLE]
where the final expectation is over the marginal distribution of under . We partition , where S_{0}:=\{j\in S:\text{j is odd}\}, S_{1+}:=\{j\in S:\text{jv_{j}=\bar{v}^{\prime}{j}}\} and S_{1-}:=\{j\in S:\text{jv{j}\neq v^{\prime}_{j}}\}. Since by construction we always have , we can apply Lemma 5 to obtain
[TABLE]
Substituting the above bound into (16), we have
[TABLE]
On the other hand, since , we also have
[TABLE]
By (17), (18) and Fano’s inequality (Yu, 1997, Lemma 3),
[TABLE]
We now choose such that \sin^{2}\gamma=\min\bigl{\{}\frac{\log|W|}{8n\lambda_{1}p\min\{1,\lambda_{1}p\}},1\bigr{\}}. Since , we obtain . Therefore,
[TABLE]
as desired. ∎
Proof of Proposition 1.
For notational simplicity, we write and for any . For , let denote the th row of . For any , we have and
[TABLE]
Thus
[TABLE]
say, where depends only on and . Note that the inequality above holds for all . Writing for convenience, we have found that . Let be the orthogonal complements of and respectively, so that and . We wish to apply Cai and Zhang (2018a, Theorem 1). To this end, note that
[TABLE]
Hence, writing , we have by Weyl’s inequality that
[TABLE]
Now, writing , we have
[TABLE]
In addition, by Cauchy–Schwarz and Jensen’s inequality,
[TABLE]
Similarly,
[TABLE]
Hence there exists , depending only on and , such that whenever , we have
[TABLE]
Now let be an SVD of . We can now apply Cai and Zhang (2018a, Theorem 1) to deduce that for such ,
[TABLE]
say. Similarly,
[TABLE]
We are now in a position to show contraction in terms of two-to-infinity norm. By Cape, Tang and Priebe (2018, Theorem 3.7),
[TABLE]
say. Note that
[TABLE]
Hence,
[TABLE]
Moreover,
[TABLE]
Finally,
[TABLE]
Write
[TABLE]
for simplicity, so that . Given that , substituting the bounds for into (7) yields that
[TABLE]
as desired. ∎
Proof of Theorem 3.
We prove this result by induction on . The case is true by definition of , so suppose that the conclusion holds for some . We make the following two claims:
- (a)
; 2. (b)
The error is further contracted by refinement, i.e., .
To prove claim (a), notice that for each , by Weyl’s inequality and the inductive hypothesis,
[TABLE]
Now, for ,
[TABLE]
On the other hand, if and , then
[TABLE]
Hence, if we choose , then , so for ,
[TABLE]
moreover, for ,
[TABLE]
Claim (a) follows. As for claim (b), note that . Taking from Proposition 1, and reducing if necessary so that , we may apply this proposition to deduce that whenever
- (i)
; 2. (ii)
,
we have . But the conditions (i) and (ii) are ensured by the inductive hypothesis and our assumptions, so the conclusion follows. ∎
Proof of Theorem 4.
Let . By Cape, Tang and Priebe (2018, Theorem 3.7), when , we have that
[TABLE]
Note that if satisfies (27) for some , then . In fact, since , there exists such that (7) implies (27), which, together with (26) ensures that .
To bound , we have
[TABLE]
where the second inequality is due to the fact that
[TABLE]
We use a covering argument to bound the supremum term. Let be a -net of the Euclidean sphere , i.e., for any , there exists a point such that . Note that for any ,
[TABLE]
which further implies that
[TABLE]
We then argue similarly as in (23), with taking the role of there (since , we correspondingly have taking the role of the incoherence parameter there), to obtain that for any ,
[TABLE]
By Vershynin (2012, Lemma 5.2), . Hence, by (7), (21) and a union bound, we have for any that
[TABLE]
Next we bound . Note that
[TABLE]
For , let and . Then
[TABLE]
By applying both parts of Lemma 1, for any and , we have . Applying Bernstein’s inequality (Boucheron, Lugosi and Massart, 2013, Theorem 2.10) yields that for any ,
[TABLE]
Therefore, a union bound with yields that
[TABLE]
Now we bound . We have that
[TABLE]
Finally, satisfies
[TABLE]
Since \|\sin\Theta(\widetilde{}\mathbf{V}_{K},\mathbf{V}_{K})\|_{\mathrm{op}}\leq\min\bigl{\{}L(\widetilde{}\mathbf{V}_{K},\mathbf{V}_{K}),1\bigr{\}}, combining our bounds for yields that there exists such that for any ,
[TABLE]
It therefore follows from Proposition 2, which applies because condition (7) for a suitable implies (8), together with the facts that and that the first conclusion of the theorem holds. The second conclusion then follows immediately. ∎
Proof of Proposition 2.
In this proof, we use the shorthand for . We represent under the orthonormal basis as follows:
[TABLE]
Define
[TABLE]
In the sequel, we regard as a corrupted version of with the off-diagonal blocks and as perturbations. We have
[TABLE]
We control the right-hand side through a concentration inequality, and for let denote the th column of . For any and ,
[TABLE]
where denotes the th column of . Note that
[TABLE]
Thus for any vector , we have by Lemma 1 that
[TABLE]
For , let . Now for any ,
[TABLE]
where the penultimate inequality uses the fact that , and the last inequality is due to Stirling’s approximation. Hence,
[TABLE]
Thus by (7) and Bernstein’s inequality (Boucheron, Lugosi and Massart, 2013, Theorem 2.10), we have that for any ,
[TABLE]
By a union bound over , for any ,
[TABLE]
Now we provide a condition under which , which ensures that is the top eigenspace of . Note that
[TABLE]
and
[TABLE]
This implies that if , then
[TABLE]
In the following, we derive an exponential tail bound for . Let and note that . Thus, for any and any integer ,
[TABLE]
Now, for each , and ,
[TABLE]
where the last inequality is due to the fact that \bigl{\|}\lVert\mathbf{y}_{i}\rVert_{\infty}\bigr{\|}_{\psi_{2}}\leq(CM\log d)^{1/2} by (10). Therefore,
[TABLE]
where . Since the above inequality holds for all , we have
[TABLE]
By a version of the Matrix Bernstein inequality for non-central absolute moments, which we give as Lemma 3, there exists a universal constant such that for any ,
[TABLE]
Now let
[TABLE]
From (25) and (26), we deduce that for any , if
[TABLE]
then \mathbb{P}^{{\boldsymbol{\Omega}}}({\cal A}^{\mathrm{c}})\leq\mathbb{P}^{{\boldsymbol{\Omega}}}\bigl{\{}\|\widetilde{}\mathbf{G}-\boldsymbol{\Sigma}_{\mathbf{y}}\|_{\mathrm{op}}\geq\lambda_{K}/4\bigr{\}}\leq 4d^{-(\xi-1)}. The desired result follows immediately by combining this with (7) and applying Yu, Wang and Samworth (2015, Theorem 2).∎
8 Auxiliary results
Lemma 1**.**
Let and be two sub-Gaussian random variables. Then we have and .
Proof.
For any , let . According to the definitions of the -norm and -norm, we have that
[TABLE]
where the penultimate inequality is due to Jensen’s inequality and the last inequality is due to the fact that . For the second inequality,
[TABLE]
as required. ∎
Lemma 2**.**
If are independent centred random variables with , then there exists a universal constant such that
[TABLE]
Proof.
Write and . From Vershynin (2012, Lemma 5.15), there exist universal constants such that for ,
[TABLE]
Setting in the above expression, the right-hand side is bounded above by . The desired result follows from the fact that (5.15) and (5.16) in Vershynin (2012) are two definitions that yield equivalent -norms. ∎
The following lemma provides a variant of the existing matrix Bernstein inequality (Tropp, 2012, Theorem 6.2). The primary difference is that we impose non-central absolute moment inequalities, as opposed to central moment inequalities. We believe that this inequality may be of independent interest, with applications beyond the scope of this paper. Recall that if , with eigendecomposition for some orthogonal , then we write .
Lemma 3** (Matrix Bernstein inequality with non-central absolute moment conditions).**
Let be independent symmetric random matrices. Assume that
[TABLE]
for some and deterministic -dimensional symmetric matrices . Define the variance parameter
[TABLE]
Then for each ,
[TABLE]
Proof.
Let be independent random matrices and variables, independent of , satisfying and for . Write and . Given , let be a leading unit-length eigenvector of . Let denote orthonormal eigenvectors of with corresponding eigenvalues ; fix , and let for . Since , we have by Jensen’s inequality that for ,
[TABLE]
We deduce that \mathbb{E}\bigl{\{}(\mathbf{v}^{\top}\widetilde{\mathbf{X}}_{i}\mathbf{v})_{+}^{q}\bigr{\}}\leq\mathbb{E}\bigl{\{}|\mathbf{v}^{\top}\widetilde{\mathbf{X}}_{i}\mathbf{v}|^{q}\bigr{\}}\leq\frac{q!}{2}R^{q-2}\mathbf{v}^{\top}\mathbf{A}_{i}^{2}\mathbf{v} for , so by Bernstein’s inequality (Boucheron, Lugosi and Massart, 2013, Corollary 2.11),
[TABLE]
We may assume that the right-hand side of the above inequality is at most , since otherwise the lemma is trivially true. Therefore,
[TABLE]
where we have used the fact that for all in the penultimate inequality.
Since and \mathbb{E}\bigl{\{}(\epsilon_{i}\mathbf{X}_{i})^{q}\bigr{\}}\preceq\mathbb{E}(|\mathbf{X}_{i}|^{q})\preceq\frac{q!}{2}R^{q-2}\mathbf{A}_{i}^{2} for , applying the matrix Bernstein inequality (Tropp, 2012, Theorem 6.2) to the sequence yields
[TABLE]
We attain the conclusion by combining the above inequality with (28). ∎
Lemma 4**.**
Let be independent random variables and let . When and , we have
[TABLE]
Proof.
By Bernstein’s inequality (van der Vaart and Wellner, 1996, Lemma 2.2.9) and a union bound,
[TABLE]
Setting , we have
[TABLE]
where we have used and in the final inequality. ∎
Lemma 5**.**
Suppose that and . Let and . Then
[TABLE]
Proof.
Since , the matrices and share the same set of eigenvalues. Hence and we have
[TABLE]
Now, by the Sherman–Morrison formula,
[TABLE]
and thus we have
[TABLE]
as required. ∎
Theorem 4 and Proposition 2 exhibit bounds on and given a deterministic observation scheme. To provide some intuition on the size of these bounds under the -homogeneous missingness setting described in Section 2.1, the following lemma derives probabilistic bounds for various norms of .
Lemma 6**.**
Assume (A5). Then there exists an event of probability at least on which each of the following inequalities hold:
- (i)
; 2. (ii)
; 3. (iii)
; 4. (iv)
; 5. (v)
.
Proof.
Define an event
[TABLE]
For , write . Then by a union bound and Bernstein’s inequality (Wainwright, 2019, Proposition 2.14), we have
[TABLE]
Note that on , we have . The desired bounds then follow respectively from the following inequalities: , , , and . ∎
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Anderson (1957) Anderson, T. W. (1957) Maximum likelihood estimates for a multivariate normal distribution when some observations are missing. J. Amer. Statist. Assoc. , 52 , 200–203.
- 2Belloni, Rosenbaum and Tsybakov (2017) Belloni, A., Rosenbaum, M. and Tsybakov, A. B. (2017) Linear and conic programming estimators in high dimensional errors-in-variables models. J. Roy. Statist. Soc., Ser. B , 79 , 939–956.
- 3Boucheron, Lugosi and Massart (2013) Boucheron, S., Lugosi, G. and Massart, P. (2013) Concentration Inequalities: A Nonasymptotic Theory of Independence . Oxford University Press, Oxford.
- 4Cai and Zhang (2018 a) Cai, T. T. and Zhang, A. (2018 a) Rate-optimal perturbation bounds for singular subspaces with applications to high-dimensional statistics. Ann. Statist. , 46 , 60–89.
- 5Cai and Zhang (2018 b) Cai, T. T. and Zhang, L. (2018 b) High-dimensional linear discriminant analysis: optimality, adaptive algorithm, and missing data. ar Xiv:1804.03018 .
- 6Candès et al. (2011) Candès, E. J., Li, X., Ma, Y. and Wright, J. (2011) Robust principal component analysis? J. ACM , 58 , 11:1–11:37.
- 7Candès and Plan (2010) Candès, E. J. and Plan, Y. (2010) Matrix completion with noise. Proc. IEEE , 98 , 925–936.
- 8Candès and Recht (2009) Candès, E. J. and Recht, B. (2009) Exact matrix completion via convex optimization. Found. Comput. Math. , 9 , 717–772.