Asymptotic performance of PCA for high-dimensional heteroscedastic data

TL;DR

This paper provides an asymptotic analysis of PCA's performance on high-dimensional heteroscedastic data, revealing that heteroscedastic noise degrades PCA recovery compared to homoscedastic noise, even with the same average noise level.

Contribution

It derives simplified asymptotic expressions for PCA performance on heteroscedastic data and demonstrates that heteroscedastic noise worsens PCA recovery relative to homoscedastic noise.

Findings

Heteroscedastic noise reduces PCA recovery performance compared to homoscedastic noise.

Simplified formulas enable easy assessment of PCA performance in high-dimensional settings.

Average noise variance overestimates PCA effectiveness for heteroscedastic data.

Abstract

Principal Component Analysis (PCA) is a classical method for reducing the dimensionality of data by projecting them onto a subspace that captures most of their variation. Effective use of PCA in modern applications requires understanding its performance for data that are both high-dimensional and heteroscedastic. This paper analyzes the statistical performance of PCA in this setting, i.e., for high-dimensional data drawn from a low-dimensional subspace and degraded by heteroscedastic noise. We provide simplified expressions for the asymptotic PCA recovery of the underlying subspace, subspace amplitudes and subspace coefficients; the expressions enable both easy and efficient calculation and reasoning about the performance of PCA. We exploit the structure of these expressions to show that, for a fixed average noise variance, the asymptotic recovery of PCA for heteroscedastic data is always worse than that for homoscedastic data (i.e., for noise variances that are equal across samples). Hence, while average noise variance is often a practically convenient measure for the overall quality of data, it gives an overly optimistic estimate of the performance of PCA for heteroscedastic data.

Equations

$$y_i=\mathbf{U\Theta }{z}_i+{\eta }_i{\epsilon }_i=\sum _{j=1}^k{\theta }_j{u}_j(z_i^{(j)}{)}^{\ast }+{\eta }_i{\epsilon }_i.$$

$${\hat{\theta }}_i^2\stackrel{a.s.}{⟶}\frac{1}{c}\max (\alpha ,{\beta }_i)\left\{1+c\sum _{\ell =1}^L\frac{p_{\ell }{\sigma }_{\ell }^2}{\max (\alpha ,{\beta }_i)-{\sigma }_{\ell }^2}\right\},$$

$$A(x)=1-c\sum _{\ell =1}^L\frac{p_{\ell }{\sigma }_{\ell }^4}{(x-{\sigma }_{\ell }^2{)}^2},B_i(x)=1-c{\theta }_i^2\sum _{\ell =1}^L\frac{p_{\ell }}{x-{\sigma }_{\ell }^2}.$$

$${\left|⟨{\hat{u}}_i,\mathrm{Span}\{u_j:{\theta }_j={\theta }_i\}⟩\right|}^2\stackrel{a.s.}{⟶}\frac{A({\beta }_i)}{{\beta }_i{B}_i^{\prime }({\beta }_i)},{\left|⟨{\hat{u}}_i,\mathrm{Span}\{u_j:{\theta }_j\ne {\theta }_i\}⟩\right|}^2\stackrel{a.s.}{⟶}0,$$

$${\left|⟨\frac{{\hat{z}}^{(i)}}{\sqrt{n}},\mathrm{Span}\{z^{(j)}:{\theta }_j={\theta }_i\}⟩\right|}^2$$

$${\left|⟨\frac{{\hat{z}}^{(i)}}{\sqrt{n}},\mathrm{Span}\{z^{(j)}:{\theta }_j\ne {\theta }_i\}⟩\right|}^2$$

$$\sum _{j:{\theta }_j={\theta }_i}⟨{\hat{u}}_i,u_j⟩{⟨\frac{{\hat{z}}^{(i)}}{\sqrt{n}},\frac{z^{(j)}}{\Vert z^{(j)}\Vert }⟩}^{\ast }\stackrel{a.s.}{⟶}\frac{A({\beta }_i)}{\sqrt{c{\beta }_i\{{\beta }_i+(1-c){\theta }_i^2\}}{B}_i^{\prime }({\beta }_i)}.$$

$$|⟨{\hat{u}}_i,\mathrm{Span}\{u_1,\dots ,u_k\}⟩{|}^2\stackrel{a.s.}{⟶}0,{\left|⟨\frac{{\hat{z}}^{(i)}}{\sqrt{n}},\mathrm{Span}\{z^{(1)},\dots ,z^{(k)}\}⟩\right|}^2\stackrel{a.s.}{⟶}0.$$

$${\hat{\theta }}_i^2\stackrel{a.s.}{⟶}\{\begin{array}{ll}{\theta }_i^2\{1+{\sigma }^2/(c{\theta }_i^2)\}(1+{\sigma }^2/{\theta }_i^2)& \text{if }c{\theta }_i^4>{\sigma }^4,\\ {\sigma }^2(1+1/\sqrt{c}{)}^2& \text{otherwise}.\end{array}$$

$${\left|⟨{\hat{u}}_i,\mathrm{Span}\{u_j:{\theta }_j={\theta }_i\}⟩\right|}^2$$

$${\left|⟨\frac{{\hat{z}}^{(i)}}{\sqrt{n}},\mathrm{Span}\{z^{(j)}:{\theta }_j={\theta }_i\}⟩\right|}^2$$

$${\left|⟨{\hat{u}}_i,\mathrm{Span}\{u_j:{\theta }_j={\theta }_i\}⟩\right|}^2$$

$${\left|⟨\frac{{\hat{z}}^{(i)}}{\sqrt{n}},\mathrm{Span}\{z^{(j)}:{\theta }_j={\theta }_i\}⟩\right|}^2$$

$$0=B_i({\beta }_i)=1-c{\theta }_i^2\sum _{\ell =1}^L\frac{p_{\ell }}{{\beta }_i-{\sigma }_{\ell }^2}<1-c{\theta }_i^2\frac{1}{{\beta }_i}.$$

$$\frac{{\hat{\theta }}_i^2}{{\theta }_i^2}\stackrel{a.s.}{⟶}\frac{{\beta }_i}{c{\theta }_i^2}\left(1+c\sum _{\ell =1}^L\frac{p_{\ell }{\sigma }_{\ell }^2}{{\beta }_i-{\sigma }_{\ell }^2}\right)$$

$$=\frac{{\beta }_i}{c{\theta }_i^2}\left(1+{\beta }_{i}c\sum _{\ell =1}^L\frac{p_{\ell }}{{\beta }_i-{\sigma }_{\ell }^2}-c\right)=\frac{{\beta }_i}{c{\theta }_i^2}\left[1+\frac{{\beta }_i}{{\theta }_i^2}\{1-B_i({\beta }_i)\}-c\right]$$

$$=\frac{{\beta }_i}{c{\theta }_i^2}\left(1+\frac{{\beta }_i}{{\theta }_i^2}-c\right)=1+\left(\frac{{\beta }_i}{c{\theta }_i^2}-1\right)\left(\frac{{\beta }_i}{{\theta }_i^2}+1\right),$$

$$\frac{1}{k}\Vert {\hat{\mathbf{U}}}^{\mathsf{H}}\mathbf{U}{\Vert }_F^2\stackrel{a.s.}{⟶}\frac{1}{k}\sum _{i=1}^k\frac{A\left({\beta }_i\right)}{{\beta }_i{B}_i^{\prime }\left({\beta }_i\right)},$$

$$\frac{1}{n}\sum _{i=1}^n$$

$$\frac{1}{k}\Vert {\hat{\mathbf{U}}}^{\mathsf{H}}\mathbf{U}{\Vert }_F^2=\frac{1}{k}\sum _{i=1}^k{\Vert {\hat{u}}_i^{\mathsf{H}}{\mathbf{U}}_{j:{\theta }_j={\theta }_i}\Vert }_2^2+{\Vert {\hat{u}}_i^{\mathsf{H}}{\mathbf{U}}_{j:{\theta }_j\ne {\theta }_i}\Vert }_2^2,$$

$$\frac{1}{n}\sum _{i=1}^n{\Vert \mathbf{U\Theta }{z}_i-\hat{\mathbf{U}}\hat{\mathbf{\Theta }}{\hat{z}}_i\Vert }_2^2$$

$$=\sum _{i=1}^k{\theta }_i^2\left[{\Vert \frac{z^{(i)}}{\sqrt{n}}\Vert }_2^2+\frac{{\hat{\theta }}_i^2}{{\theta }_i^2}-2\Re \left\{\frac{{\hat{\theta }}_i}{{\theta }_i}\sum _{j=1}^k\frac{{\theta }_j}{{\theta }_i}⟨{\hat{u}}_i,u_j⟩{⟨\frac{{\hat{z}}^{(i)}}{\sqrt{n}},\frac{z^{(j)}}{\sqrt{n}}⟩}^{\ast }\right\}\right],$$

$$\sum _{j:{\theta }_j={\theta }_i}\frac{{\theta }_j}{{\theta }_i}⟨{\hat{u}}_i,u_j⟩{⟨\frac{{\hat{z}}^{(i)}}{\sqrt{n}},\frac{z^{(j)}}{\sqrt{n}}⟩}^{\ast }=\sum _{j:{\theta }_j={\theta }_i}(1)⟨{\hat{u}}_i,u_j⟩{⟨\frac{{\hat{z}}^{(i)}}{\sqrt{n}},\frac{z^{(j)}}{\sqrt{n}}⟩}^{\ast }$$

$${\hat{\theta }}_i^2\stackrel{a.s.}{⟶}{\theta }_i^2+{\theta }_i^2\left(\frac{{\beta }_i}{c{\theta }_i^2}-1\right)\left(\frac{{\beta }_i}{{\theta }_i^2}+1\right)\ge {\theta }_i^2\left(1+\frac{{\bar{\sigma }}^2}{c{\theta }_i^2}\right)\left(1+\frac{{\bar{\sigma }}^2}{{\theta }_i^2}\right),$$

$${\left|⟨{\hat{u}}_i,\mathrm{Span}\{u_j:{\theta }_j={\theta }_i\}⟩\right|}^2\stackrel{a.s.}{⟶}\frac{A\left({\beta }_i\right)}{{\beta }_i{B}_i^{\prime }\left({\beta }_i\right)}\le \frac{c-{\bar{\sigma }}^4/{\theta }_i^4}{c+{\bar{\sigma }}^2/{\theta }_i^2},$$

$${\left|⟨\frac{{\hat{z}}^{(i)}}{\sqrt{n}},\mathrm{Span}\{z^{(j)}:{\theta }_j={\theta }_i\}⟩\right|}^2\stackrel{a.s.}{⟶}\frac{A({\beta }_i)}{c\{{\beta }_i+(1-c){\theta }_i^2\}{B}_i^{\prime }({\beta }_i)}\le \frac{c-{\bar{\sigma }}^4/{\theta }_i^4}{c(1+{\bar{\sigma }}^2/{\theta }_i^2)},$$

$$\frac{1}{c{\theta }_i^2}=\sum _{\ell =1}^L\frac{p_{\ell }}{x-{\sigma }_{\ell }^2}.$$

$$\mathbf{Y}=(y_1,\dots ,y_n)=\mathbf{U\Theta }{\mathbf{Z}}^{\mathsf{H}}+\mathbf{EH}\in {\mathbb{C}}^{d\times n},$$

$$\mathbf{R}=[\begin{array}[]{cc}\mathbf{\tilde{U}}&\mathbf{\tilde{U}}^{\perp}\end{array}][\begin{array}[]{cc}\mathbf{U}&\mathbf{U}^{\perp}\end{array}]^{\mathsf{H}}=\mathbf{\tilde{U}}\mathbf{U}^{\mathsf{H}}+\mathbf{\tilde{U}}^{\perp}(\mathbf{U}^{\perp})^{\mathsf{H}},$$

$$\tilde{\mathbf{Y}}=\frac{1}{\sqrt{n}}\mathbf{RY}.$$

$$\tilde{\mathbf{Y}}=\mathbf{P}+\mathbf{X},$$

$$\mathbf{P}=\frac{1}{\sqrt{n}}\mathbf{R}\left(\mathbf{U\Theta }{\mathbf{Z}}^{\mathsf{H}}\right)=\frac{1}{\sqrt{n}}\tilde{\mathbf{U}}\mathbf{\Theta }{\mathbf{Z}}^{\mathsf{H}}=\sum _{i=1}^k{\theta }_i{\tilde{u}}_i{\left(\frac{1}{\sqrt{n}}{z}^{\left(i\right)}\right)}^{\mathsf{H}},\mathbf{X}=\frac{1}{\sqrt{n}}\mathbf{R}\left(\mathbf{EH}\right)=\left(\frac{1}{\sqrt{n}}\mathbf{RE}\right)\mathbf{H}.$$

$$|⟨\mathbf{R}{\hat{u}}_i,\mathrm{Span}\{{\tilde{u}}_j:{\theta }_j={\theta }_i\}⟩{|}^2$$

Full text

††journal: Journal of Multivariate Analysis

1 Introduction

Principal Component Analysis (PCA) is a classical method for reducing the dimensionality of data by representing them in terms of a new set of variables, called principal components, where variation in the data is largely captured by the first few principal components [23]. This paper analyzes the asymptotic performance of PCA for data with heteroscedastic noise. In particular, we consider the classical and commonly employed unweighted form of PCA that treats all samples equally and remains a natural choice in applications where estimates of the noise variances are unavailable or one hopes the noise is “close enough” to being homoscedastic. Our analysis uncovers several practical new insights for this setting; the findings both broaden our understanding of PCA and also precisely characterize the impact of heteroscedasticity.

Given zero-mean sample vectors $y_{1},\dots,y_{n}\in\mathbb{C}^{d}$, the first $k$ principal components $\hat{u}_{1},\dots,\hat{u}_{k}\in\mathbb{C}^{d}$ and corresponding squared PCA amplitudes $\hat{\theta}_{1}^{2},\dots,\hat{\theta}_{k}^{2}\in\mathbb{R}_{+}$ are the first $k$ eigenvectors and eigenvalues, respectively, of the sample covariance matrix $(y_{1}y_{1}^{\mathsf{H}}+\cdots+y_{n}y_{n}^{\mathsf{H}})/n$. The associated score vectors $\hat{z}^{(1)},\dots,\hat{z}^{(k)}\in\mathbb{C}^{n}$ are standardized projections given, for each $i\in\{1,\ldots,k\}$, by $\hat{z}^{(i)}=(1/\hat{\theta}_{i})\{\hat{u}_{i}^{\mathsf{H}}(y_{1},\ldots,y_{n})\}^{\mathsf{H}}$. The principal components $\hat{u}_{1},\ldots,\hat{u}_{k}$, PCA amplitudes $\hat{\theta}_{1},\ldots,\hat{\theta}_{k}$ and score vectors $\hat{z}^{(1)},\dots,\hat{z}^{(k)}$ are efficiently obtained from the data matrix $(y_{1},\ldots,y_{n})\in\mathbb{C}^{d\times n}$ as its left singular vectors, (scaled) singular values and (scaled) right singular vectors, respectively.

A natural setting for PCA is when data are noisy measurements of points drawn from a subspace. In this case, the first few principal components $\hat{u}_{1},\ldots,\hat{u}_{k}$ form an estimated basis for the underlying subspace; if they recover the underlying subspace accurately then the low-dimensional scores $\hat{z}^{(1)},\ldots,\hat{z}^{(k)}$ will largely capture the meaningful variation in the data. This paper analyzes how well the first $k$ principal components $\hat{u}_{1},\ldots,\hat{u}_{k}$, PCA amplitudes $\hat{\theta}_{1},\ldots,\hat{\theta}_{k}$ and score vectors $\hat{z}^{(1)},\ldots,\hat{z}^{(k)}$ recover their underlying counterparts when the data are heteroscedastic, that is, when the noise in the data has non-uniform variance across samples.

1.1 High-dimensional, heteroscedastic data

Dimensionality reduction is a fundamental task, so PCA has been applied in a broad variety of both traditional and modern settings. See [23] for a thorough review of PCA and some of its important traditional applications. A sample of modern application areas include medical imaging [2, 33], classification for cancer data [35], anomaly detection on computer networks [24], environmental monitoring [31, 39] and genetics [25], to name just a few.

It is common in modern applications in particular for the data to be high-dimensional (i.e., the number of variables measured is comparable with or larger than the number of samples), which has motivated the development of new techniques and theory for this regime [22]. It is also common for modern data sets to have heteroscedastic noise. For example, Cochran and Horne [11] apply a PCA variant to spectrophotometric data from the study of chemical reaction kinetics. The spectrophotometric data are absorptions at various wavelengths over time, and measurements are averaged over increasing windows of time causing the amount of noise to vary across time. Another example is given in [36], where data are astronomical measurements of stars taken at various times; here changing atmospheric effects cause the amount of noise to vary across time. More generally, in the era of big data where inference is made using numerous data samples drawn from a myriad of different sources, one can expect that both high-dimensionality and heteroscedasticity will be the norm. It is important to understand the performance of PCA in such settings.

1.2 Contributions of this paper

This paper provides simplified expressions for the performance of PCA from heteroscedastic data in the limit as both the number of samples and dimension tend to infinity. The expressions quantify the asymptotic recovery of an underlying subspace, subspace amplitudes and coefficients by the principal components, PCA amplitudes and scores, respectively. The asymptotic recoveries are functions of the samples per ambient dimension, the underlying subspace amplitudes and the distribution of noise variances. Forming the expressions involves first connecting several results from random matrix theory [3, 5] to obtain initial expressions for asymptotic recovery that are difficult to evaluate and analyze, and then exploiting a nontrivial structure in the expressions to obtain much simpler algebraic descriptions. These descriptions enable both easy and efficient calculation and reasoning about the asymptotic performance of PCA.

The impact of heteroscedastic noise, in particular, is not immediately obvious given results of prior literature. How much do a few noisy samples degrade the performance of PCA? Is heteroscedasticity ever beneficial for PCA? Our simplified expressions enable such questions to be answered. In particular, we use these expressions to show that, for a fixed average noise variance, the asymptotic subspace recovery, amplitude recovery and coefficient recovery are all worse for heteroscedastic data than for homoscedastic data (i.e., for noise variances that are equal across samples), confirming a conjecture in [18]. Hence, while average noise variance is often a practically convenient measure for the overall quality of data, it gives an overly optimistic estimate of PCA performance. This analysis provides a deeper understanding of how PCA performs in the presence of heteroscedastic noise.

1.3 Relationship to previous works

Homoscedastic noise has been well-studied, and there are many nice results characterizing PCA in this setting. Benaych-Georges and Nadakuditi [5] give an expression for asymptotic subspace recovery, also found in [21, 29, 32], in the limit as both the number of samples and ambient dimension tend to infinity. As argued in [21], the expression in [5] reveals that asymptotic subspace recovery is perfect only when the number of samples per ambient dimension tends to infinity, so PCA is not (asymptotically) consistent for high-dimensional data. Various alternatives [6, 15, 21] can regain consistency by exploiting sparsity in the covariance matrix or in the principal components. As discussed in [5, 29], the expression in [5] also exhibits a phase transition: the number of samples per ambient dimension must be sufficiently high to obtain non-zero subspace recovery (i.e., for any subspace recovery to occur). This paper generalizes the expression in [5] to heteroscedastic noise; homoscedastic noise is a special case and is discussed in Section 2.3. Once again, (asymptotic) consistency is obtained when the number of samples per ambient dimension tends to infinity, and there is a phase transition between zero recovery and non-zero recovery.

PCA is known to generally perform well in the presence of low to moderate homoscedastic noise and in the presence of missing data [10]. When the noise is standard normal, PCA gives the maximum likelihood (ML) estimate of the subspace [37]. In general, [37] proposes finding the ML estimate via expectation maximization. Conventional PCA is not an ML estimate of the subspace for heteroscedastic data, but it remains a natural choice in applications where we might expect noise to be heteroscedastic but hope it is “close enough” to being homoscedastic. Even for mostly homoscedastic data, however, PCA performs poorly when the heteroscedasticity is due to gross errors (i.e., outliers) [13, 19, 23], which has motivated the development and analysis of robust variants; see [8, 9, 12, 16, 17, 26, 34, 40, 42] and their corresponding bibliographies. This paper provides expressions for asymptotic recovery that enable rigorous understanding of the impact of heteroscedasticity.

The generalized spiked covariance model, proposed and analyzed in [4] and [41], generalizes homoscedastic noise in an alternate way. It extends the Johnstone spiked covariance model [20, 21] (a particular homoscedastic setting) by using a population covariance that allows, among other things, non-uniform noise variances within each sample. Non-uniform noise variances within each sample may arise, for example, in applications where sample vectors are formed by concatenating the measurements of intrinsically different quantities. This paper considers data with non-uniform noise variances across samples instead; we model noise variances within each sample as uniform. Data with non-uniform noise variances across samples arise, for example, in applications where samples come from heterogeneous sources, some of which are better quality (i.e., lower noise) than others. See Section S1 of the Online Supplement for a more detailed discussion of connections to spiked covariance models.

Our previous work [18] analyzed the subspace recovery of PCA for heteroscedastic noise but was limited to real-valued data coming from a random one-dimensional subspace where the number of samples exceeded the data dimension. This paper extends that analysis to the more general setting of real- or complex-valued data coming from a deterministic low-dimensional subspace where the number of samples no longer needs to exceed the data dimension. This paper also extends the analysis of [18] to include the recovery of the underlying subspace amplitudes and coefficients. In both works, we use the main results of [5] to obtain initial expressions relating asymptotic recovery to the limiting noise singular value distribution.

The main results of [29] provide non-asymptotic results (i.e., probabilistic approximation results for finite samples in finite dimension) for homoscedastic noise limited to the special case of one-dimensional subspaces. Signal-dependent noise was recently considered in [38], where they analyze the performance of PCA and propose a new generalization of PCA that performs better in certain regimes. A recent extension of [5] to linearly reduced data is presented in [14] and may be useful for analyzing weighted variants of PCA. Such analyses are beyond the scope of this paper, but are interesting avenues for further study.

1.4 Organization of the paper

Section 2 describes the model we consider and states the main results: simplified expressions for asymptotic PCA recovery and the fact that PCA performance is best (for a fixed average noise variance) when the noise is homoscedastic. Section 3 uses the main results to provide a qualitative analysis of how the model parameters (e.g., samples per ambient dimension and the distribution of noise variances) affect PCA performance under heteroscedastic noise. Section 4 compares the asymptotic recovery with non-asymptotic (i.e., finite) numerical simulations. The simulations demonstrate good agreement as the ambient dimension and number of samples grow large; when these values are small the asymptotic recovery and simulation differ but have the same general behavior. Sections 5 and 6 prove the main results. Finally, Section 7 discusses the findings and describes avenues for future work.

2 Main results

2.1 Model for heteroscedastic data

We model $n$ heteroscedastic sample vectors $y_{1},\ldots,y_{n}\in\mathbb{C}^{d}$ from a $k$-dimensional subspace as

[TABLE]

The following are deterministic:

• $\mathbf{U}=(u_{1},\dots,u_{k})\in\mathbb{C}^{d\times k}$ forms an orthonormal basis for the subspace,
• $\mathbf{\Theta}=\operatorname*{diag}(\theta_{1},\dots,\theta_{k})\in\mathbb{R}_{+}^{k\times k}$ is a diagonal matrix of amplitudes,
• $\eta_{i}\in\{\sigma_{1},\dots,\sigma_{L}\}$ are each one of $L$ noise standard deviations $\sigma_{1},\dots,\sigma_{L}$,

and we define $n_{1}$ to be the number of samples with $\eta_{i}=\sigma_{1}$, $n_{2}$ to be the number of samples with $\eta_{i}=\sigma_{2}$ and so on, where $n_{1}+\cdots+n_{L}=n$.

The following are random and independent:

• $z_{i}\in\mathbb{C}^{k}$ are iid sample coefficient vectors that have iid entries with mean ${\rm E}(z_{ij})=0$, variance ${\rm E}|z_{ij}|^{2}=1$, and a distribution satisfying the log-Sobolev inequality [1],
• $\varepsilon_{i}\in\mathbb{C}^{d}$ are unitarily invariant iid noise vectors that have iid entries with mean ${\rm E}(\varepsilon_{ij})=0$, variance ${\rm E}|\varepsilon_{ij}|^{2}=1$ and bounded fourth moment ${\rm E}|\varepsilon_{ij}|^{4}<\infty$,

and we define the $k$ (component) coefficient vectors $z^{(1)},\dots,z^{(k)}\in\mathbb{C}^{n}$ such that the $i$th entry of $z^{(j)}$ is $z^{(j)}_{i}=(z_{ij})^{*}$, the complex conjugate of the $j$th entry of $z_{i}$. Defining the coefficient vectors in this way is convenient for stating and proving the results that follow, as they more naturally correspond to right singular vectors of the data matrix formed by concatenating $y_{1},\dots,y_{n}$ as columns.

The model extends the Johnstone spiked covariance model [20, 21] by incorporating heteroscedasticity (see Section S1 of the Online Supplement for a detailed discussion). We also allow complex-valued data, as it is of interest in important signal processing applications such as medical imaging; for example, data obtained in magnetic resonance imaging are complex-valued.

Remark 1**.**

By unitarily invariant, we mean that left multiplication of $\varepsilon_{i}$ by any unitary matrix does not change the joint distribution of its entries. As in our previous work [18], this assumption can be dropped if instead the subspace $\mathbf{U}$ is randomly drawn according to either the “orthonormalized model” or “iid model” of [5]. Under these models, the subspace $\mathbf{U}$ is randomly chosen in an isotropic manner.

Remark 2**.**

The above conditions are satisfied, for example, when the entries $z_{ij}$ and $\varepsilon_{ij}$ are circularly symmetric complex normal $\mathcal{CN}(0,1)$ or real-valued normal $\mathcal{N}(0,1)$. Rademacher random variables (i.e., $\pm 1$ with equal probability) are another choice for coefficient entries $z_{ij}$; see Section 2.3.2 of [1] for discussion of the log-Sobolev inequality. We are unaware of non-Gaussian distributions satisfying all conditions for noise entries $\varepsilon_{ij}$, but as noted in Remark 1, unitary invariance can be dropped if we assume the subspace is randomly drawn as in [5].

Remark 3**.**

The assumption that noise entries $\varepsilon_{ij}$ are identically distributed with bounded fourth moment can be relaxed when they are real-valued as long as an aggregate of their tails still decays sufficiently quickly, i.e., as long as they satisfy Condition 1.3 from [30]. In this setting, the results of [30] replace those of [3] in the proof.

2.2 Simplified expressions for asymptotic recovery

The following theorem describes how well the PCA estimates $\hat{u}_{1},\ldots,\hat{u}_{k}$, $\hat{\theta}_{1},\ldots,\hat{\theta}_{k}$ and $\hat{z}^{(1)},\ldots,\hat{z}^{(k)}$ recover the underlying subspace basis $u_{1},\ldots,u_{k}$, subspace amplitudes $\theta_{1},\ldots,\theta_{k}$ and coefficient vectors $z^{(1)},\ldots,z^{(k)}$, as a function of the sample-to-dimension ratio $n/d\to c>0$, the subspace amplitudes $\theta_{1},\ldots,\theta_{k}$, the noise variances $\sigma_{1}^{2},\ldots,\sigma_{L}^{2}$ and corresponding proportions $n_{\ell}/n\to p_{\ell}$ for each $\ell\in\{1,\ldots,L\}$. One may generally expect performance to improve with increasing sample-to-dimension ratio and subspace amplitudes; Theorem 1 provides the precise dependence on these parameters as well as on the noise variances and their proportions.

Theorem 1** (Recovery of individual components).**

Suppose that the sample-to-dimension ratio $n/d\to c>0$ and the noise variance proportions $n_{\ell}/n\to p_{\ell}$ for $\ell\in\{1,\dots,L\}$ as $n,d\to\infty$. Then the $i$th PCA amplitude $\hat{\theta}_{i}$ is such that

[TABLE]

where $\alpha$ and $\beta_{i}$ are, respectively, the largest real roots of

[TABLE]

Furthermore, if $A(\beta_{i})>0$, then the $i$th principal component $\hat{u}_{i}$ is such that

[TABLE]

the normalized score vector $\hat{z}^{(i)}/\sqrt{n}$ is such that

[TABLE]

and

[TABLE]

Section 5 presents the proof of Theorem 1. The expressions can be easily and efficiently computed. The hardest part is finding the largest roots of the univariate rational functions $A(x)$ and $B_{i}(x)$, but off-the-shelf solvers can do this efficiently. See [18] for an example of similar calculations.

The projection $|\langle\hat{u}_{i},\operatorname*{Span}\{u_{j}:\theta_{j}=\theta_{i}\}\rangle|^{2}$ in Theorem 1 is the square cosine principal angle between the $i$th principal component $\hat{u}_{i}$ and the span of the basis elements with subspace amplitudes equal to $\theta_{i}$. When the subspace amplitudes are distinct, $|\langle\hat{u}_{i},\operatorname*{Span}\{u_{j}:\theta_{j}=\theta_{i}\}\rangle|^{2}=\left|\left\langle\hat{u}_{i},u_{i}\right\rangle\right|^{2}$ is the square cosine angle between $\hat{u}_{i}$ and $u_{i}$. This value is related by a constant to the squared error between the two (unit norm) vectors and is one among several natural performance metrics for subspace estimation. Similar observations hold for $|\langle\hat{z}^{(i)}/\sqrt{n},\operatorname*{Span}\{z^{(j)}:\theta_{j}=\theta_{i}\}\rangle|^{2}$. Note that $\hat{z}^{(i)}/\sqrt{n}$ has unit norm.

The expressions (4), (5) and (6) apply only if $A(\beta_{i})>0$. The following conjecture predicts a phase transition at $A(\beta_{i})=0$ so that asymptotic recovery is zero for $A(\beta_{i})\leq 0$.

Conjecture 1** (Phase transition).**

Suppose (as in Theorem 1) that the sample-to-dimension ratio $n/d\to c>0$ and the noise variance proportions $n_{\ell}/n\to p_{\ell}$ for $\ell\in\{1,\dots,L\}$ as $n,d\to\infty$. If $A(\beta_{i})\leq 0$, then the $i$th principal component $\hat{u}_{i}$ and the normalized score vector $\hat{z}^{(i)}/\sqrt{n}$ are such that

[TABLE]

This conjecture is true for a data model having Gaussian coefficients and homoscedastic Gaussian noise as shown in [32]. It is also true for a one-dimensional subspace (i.e., $k=1$) as we showed in [18]. Proving it in general would involve showing that the singular values of the matrix whose columns are the noise vectors exhibit repulsion behavior; see Remark 2.13 of [5].

2.3 Homoscedastic noise as a special case

For homoscedastic noise with variance $\sigma^{2}$, $A(x)=1-c\sigma^{4}/(x-\sigma^{2})^{2}$ and $B_{i}(x)=1-c\theta_{i}^{2}/(x-\sigma^{2})$. The largest real roots of these functions are, respectively, $\alpha=(1+\sqrt{c})\sigma^{2}$ and $\beta_{i}=\sigma^{2}+c\theta_{i}^{2}$. Thus the asymptotic PCA amplitude (2) becomes

[TABLE]

Further, if $c\theta_{i}^{4}>\sigma^{4}$, then the non-zero portions of asymptotic subspace recovery (4) and coefficient recovery (5) simplify to

[TABLE]

These limits agree with the homoscedastic results in [5, 7, 21, 29, 32]. As noted in Section 2.2, Conjecture 1 is known to be true when the coefficients are Gaussian and the noise is both homoscedastic and Gaussian, in which case (8) becomes

[TABLE]

See Section 2 of [21] and Section 2.3 of [32] for a discussion of this result.

2.4 Bias of the PCA amplitudes

The simplified expression in (2) enables us to immediately make two observations about the recovery of the subspace amplitudes $\theta_{1},\dots,\theta_{k}$ by the PCA amplitudes $\hat{\theta}_{1},\dots,\hat{\theta}_{k}$.

Remark 4** (Positive bias in PCA amplitudes).**

The largest real root $\beta_{i}$ of $B_{i}(x)$ is greater than $\max_{\ell}(\sigma_{\ell}^{2})$. Thus $1/(\beta_{i}-\sigma_{\ell}^{2})>1/\beta_{i}$ for $\ell\in\{1,\dots,L\}$ and so evaluating (3) at $\beta_{i}$ yields

[TABLE]

As a result, $\beta_{i}>c\theta_{i}^{2}$, so the asymptotic PCA amplitude (2) exceeds the subspace amplitude, i.e., $\hat{\theta}_{i}$ is positively biased and is thus an inconsistent estimate of $\theta_{i}$. This is a general phenomenon for noisy data and motivates asymptotically optimal shrinkage in [27].

Remark 5** (Alternate formula for amplitude bias).**

If $A(\beta_{i})\geq 0$, then $\beta_{i}\geq\alpha$ because $A(x)$ and $B_{i}(x)$ are both increasing functions for $x>\max_{\ell}(\sigma_{\ell}^{2})$. Thus, the asymptotic amplitude bias is

[TABLE]

where we have applied (2), divided the summand with respect to $\sigma_{\ell}^{2}$, used the facts that $p_{1}+\cdots+p_{L}=1$ and $B_{i}(\beta_{i})=0$, and finally factored. The expression (5) shows that the positive bias is an increasing function of $\beta_{i}$ when $A(\beta_{i})\geq 0$.

2.5 Overall subspace and signal recovery

Overall subspace recovery is more useful than individual component recovery when subspace amplitudes are equal and so individual basis elements are not identifiable. It is also more relevant when we are most interested in recovering or denoising low-dimensional signals in a subspace. Overall recovery of the low-dimensional signal, quantified here by mean square error, is useful for understanding how well PCA “denoises” the data taken as a whole.

Corollary 1** (Overall recovery).**

Suppose (as in Theorem 1) that the sample-to-dimension ratio $n/d\to c>0$ and the noise variance proportions $n_{\ell}/n\to p_{\ell}$ for $\ell\in\{1,\dots,L\}$ as $n,d\to\infty$. If $A(\beta_{1}),\dots,A(\beta_{k})>0$, then the subspace estimate $\mathbf{\hat{U}}=(\hat{u}_{1},\dots,\hat{u}_{k})\in\mathbb{C}^{d\times k}$ from PCA is such that

[TABLE]

and the mean square error is

[TABLE]

where $A(x)$, $B_{i}(x)$ and $\beta_{i}$ are as in Theorem 1, and $\hat{z}_{i}$ is the vector of score entries for the $i$th sample.

Proof of Corollary 1.

The subspace recovery can be decomposed as

[TABLE]

where the columns of $\mathbf{U}_{j:\theta_{j}=\theta_{i}}$ are the basis elements $u_{j}$ with subspace amplitude $\theta_{j}$ equal to $\theta_{i}$, and the remaining basis elements are the columns of $\mathbf{U}_{j:\theta_{j}\neq\theta_{i}}$. Asymptotic overall subspace recovery (10) follows by noting that these terms are exactly the square cosine principal angles in (4) of Theorem 1.

The mean square error can also be decomposed as

[TABLE]

where $\mathbf{Z}=(z^{(1)},\ldots,z^{(k)})\in\mathbb{C}^{n\times k}$, $\mathbf{\hat{Z}}=(\hat{z}^{(1)},\ldots,\hat{z}^{(k)})\in\mathbb{C}^{n\times k}$ and $\Re$ denotes the real part of its argument. The first term of (12) has almost sure limit $1$ by the law of large numbers. The almost sure limit of the second term is obtained from (5). We can disregard the summands in the inner sum for which $\theta_{j}\neq\theta_{i}$; by (4) and (5) these terms have an almost sure limit of zero (the inner products both vanish). The rest of the inner sum

[TABLE]

has the same almost sure limit as in (6) because $\|z^{(i)}/\sqrt{n}\|^{2}\to 1$ as $n\to\infty$. Combining these almost sure limits and simplifying yields (11). ∎

2.6 Importance of homoscedasticity

How important is homoscedasticity for PCA? Does having some low noise data outweigh the cost of introducing heteroscedasticity? Consider the following three settings:

• 1.-

All samples have noise variance $1$ (i.e., data are homoscedastic).

• 2.-

$99\%$ of samples have noise variance $1.01$ but $1\%$ have noise variance $0.01$.

• 3.-

$99\%$ of samples have noise variance $0.01$ but $1\%$ have noise variance $99.01$.

In all three settings, the average noise variance is $1$. We might expect PCA to perform well in Setting 1 because it has the smallest maximum noise variance. However, Setting 2 may seem favorable because we obtain samples with very small noise, and suffer only a slight increase in noise for the rest. Setting 3 may seem favorable because most of the samples have very small noise. However, we might also expect PCA to perform poorly because $1\%$ of samples have very large noise and will likely produce gross errors (i.e., outliers). Between all three, it is not initially obvious what setting PCA will perform best in. The following theorem shows that PCA performs best when the noise is homoscedastic, as in Setting 1.

Theorem 2**.**

Homoscedastic noise produces the best asymptotic PCA amplitude (2), subspace recovery (4) and coefficient recovery (5) in Theorem 1 for a given average noise variance $\bar{\sigma}^{2}=p_{1}\sigma_{1}^{2}+\cdots+p_{L}\sigma_{L}^{2}$ over all distributions of noise variances for which $A(\beta_{i})>0$. Namely, homoscedastic noise minimizes (2) (and hence the positive bias) and it maximizes (4) and (5).

Concretely, suppose we had $c=10$ samples per dimension and a subspace amplitude of $\theta_{i}=1$. Then the asymptotic subspace recoveries (4) given in Theorem 1 evaluate to $0.818$ in Setting 1, $0.817$ in Setting 2 and [math] in Setting 3; asymptotic recovery is best in Setting 1 as predicted by Theorem 2. Recovery is entirely lost in Setting 3, consistent with the observation that PCA is not robust to gross errors. In Setting 2, only using the $1\%$ of samples with noise variance $0.01$ (resulting in $0.1$ samples per dimension) yields an asymptotic subspace recovery of $0.908$ and so we may hope that recovery with all data could be better. Theorem 2 rigorously shows that PCA does not fully exploit these high quality samples and instead performs worse in Setting 2 than in Setting 1, if only slightly.

Section 6 presents the proof of Theorem 2. It is notable that Theorem 2 holds for all proportions $p$, sample-to-dimension ratios $c$ and subspace amplitudes $\theta_{i}$; there are no settings where PCA benefits from heteroscedastic noise over homoscedastic noise with the same average variance. The following corollary is equivalent and provides an alternate way of viewing the result.

Corollary 2** (Bounds on asymptotic recovery).**

If $A(\beta_{i})\geq 0$ then the asymptotic PCA amplitude (2) is bounded as

[TABLE]

the asymptotic subspace recovery (4) is bounded as

[TABLE]

and the asymptotic coefficient recovery (5) is bounded as

[TABLE]

where $\bar{\sigma}^{2}=p_{1}\sigma_{1}^{2}+\cdots+p_{L}\sigma_{L}^{2}$ is the average noise variance and the bounds are met with equality if and only if $\sigma_{1}^{2}=\cdots=\sigma_{L}^{2}$.

Proof of Corollary 2.

The bounds (13), (14) and (15) follow immediately from Theorem 2 and the expressions for homoscedastic noise (7) and (8) in Section 2.3. ∎

Corollary 2 highlights that while average noise variance may be a practically convenient measure for the overall quality of data, it can lead to an overly optimistic estimate of the performance of PCA for heteroscedastic data. The expressions (2), (4) and (5) in Theorem 1 are more accurate.

Remark 6** (Average inverse noise variance).**

Average inverse noise variance $\mathcal{I}=p_{1}\times 1/\sigma_{1}^{2}+\cdots+p_{L}\times 1/\sigma_{L}^{2}$ is another natural measure for the overall quality of data. In particular, it is the (scaled) Fisher information for heteroscedastic Gaussian measurements of a fixed scalar. Theorem 2 implies that homoscedastic noise also produces the best asymptotic PCA performance for a given average inverse noise variance; note that homoscedastic noise minimizes the average noise variance in this case. Thus, average inverse noise variance can also lead to an overly optimistic estimate of the performance of PCA for heteroscedastic data.

3 Impact of parameters

The simplified expressions in Theorem 1 for the asymptotic performance of PCA provide insight into the impact of the model parameters: sample-to-dimension ratio $c$, subspace amplitudes $\theta_{1},\dots,\theta_{k}$, proportions $p_{1},\dots,p_{L}$ and noise variances $\sigma_{1}^{2},\dots,\sigma_{L}^{2}$. For brevity, we focus on the asymptotic subspace recovery (4) of the $i$th component; similar phenomena occur for the asymptotic PCA amplitudes (2) and coefficient recovery (5) as we show in Section S3 of the Online Supplement.

3.1 Impact of sample-to-dimension ratio $c$ and subspace amplitude $\theta_{i}$

Suppose first that there is only one noise variance fixed at $\sigma_{1}^{2}=1$ while we vary the sample-to-dimension ratio $c$ and subspace amplitude $\theta_{i}$. This is the homoscedastic setting described in Section 2.3. Figure 1(a) illustrates the expected behavior: decreasing the subspace amplitude $\theta_{i}$ degrades asymptotic subspace recovery (4) but the lost performance could be regained by increasing the number of samples. Figure 1(a) also illustrates a phase transition: a sufficient number of samples with a sufficiently large subspace amplitude is necessary to have an asymptotic recovery greater than zero. Note that in all such figures, we label the axis $\left|\langle\hat{u}_{i},u_{i}\rangle\right|^{2}$ to indicate the asymptotic recovery on the right hand side of (4).

Now suppose that there are two noise variances $\sigma_{1}^{2}=0.8$ and $\sigma_{2}^{2}=1.8$ occurring in proportions $p_{1}=80\%$ and $p_{2}=20\%$. The average noise variance is still $1$, and Figure 1(b) illustrates similar overall features to the homoscedastic case. Decreasing subspace amplitude $\theta_{i}$ once again degrades asymptotic subspace recovery (4) and the lost performance could be regained by increasing the number of samples. However, the phase transition is further up and to the right compared to the homoscedastic case. This is consistent with Theorem 2; PCA performs worse on heteroscedastic data than it does on homoscedastic data of the same average noise variance, and thus more samples or a larger subspace amplitude are needed to recover the subspace basis element.

3.2 Impact of proportions $p_{1},\dots,p_{L}$

Suppose that there are two noise variances $\sigma_{1}^{2}=0.1$ and $\sigma_{2}^{2}=3.25$ occurring in proportions $p_{1}=1-p_{2}$ and $p_{2}$, where the sample-to-dimension ratio is $c=10$ and the subspace amplitude is $\theta_{i}=1$. Figure 2 shows the asymptotic subspace recovery (4) as a function of the proportion $p_{2}$. Since $\sigma_{2}^{2}$ is significantly larger, it is natural to think of $p_{2}$ as a fraction of contaminated samples. As expected, performance generally degrades as $p_{2}$ increases and low noise samples with noise variance $\sigma_{1}^{2}$ are traded for high noise samples with noise variance $\sigma_{2}^{2}$. The performance is best when $p_{2}=0$ and all the samples have the smaller noise variance $\sigma_{1}^{2}$ (i.e., there is no contamination).

It is interesting that the asymptotic subspace recovery in Figure 2 has a steeper slope initially for $p_{2}$ close to zero and then a shallower slope for $p_{2}$ close to one. Thus the benefit of reducing the contamination fraction varies across the range.

3.3 Impact of noise variances $\sigma_{1}^{2},\dots,\sigma_{L}^{2}$

Suppose that there are two noise variances $\sigma_{1}^{2}$ and $\sigma_{2}^{2}$ occurring in proportions $p_{1}=70\%$ and $p_{2}=30\%$, where the sample-to-dimension ratio is $c=10$ and the subspace amplitude is $\theta_{i}=1$. Figure 3 shows the asymptotic subspace recovery (4) as a function of the noise variances $\sigma_{1}^{2}$ and $\sigma_{2}^{2}$. As expected, performance typically degrades with increasing noise variances. However, there is a curious regime around $\sigma_{1}^{2}=0$ and $\sigma_{2}^{2}=4$ where increasing $\sigma_{1}^{2}$ slightly from zero improves asymptotic performance; the contour lines point slightly up and to the right. We have also observed this phenomenon in finite-dimensional simulations, so this effect is not simply an asymptotic artifact. This surprising phenomenon is an interesting avenue for future exploration.

The contours in Figure 3 are generally horizontal for small $\sigma_{1}^{2}$ and vertical for small $\sigma_{2}^{2}$. This indicates that when the gap between the two largest noise variances is “sufficiently” wide, the asymptotic subspace recovery (4) is roughly determined by the largest noise variance. While initially unexpected, this property can be intuitively understood by recalling that $\beta_{i}$ is the largest value of $x$ satisfying

[TABLE]

When the gap between the two largest noise variances is wide, the largest noise variance is significantly larger than the rest and it dominates the sum in (16) for $x>\max_{\ell}(\sigma_{\ell}^{2})$, i.e., where $\beta_{i}$ occurs. Thus $\beta_{i}$, and similarly, $A(\beta_{i})$ and $B_{i}^{\prime}(\beta_{i})$ are roughly determined by the largest noise variance.

The precise relative impact of each noise variance $\sigma_{\ell}^{2}$ depends on its corresponding proportion $p_{\ell}$, as shown by the asymmetry of Figure 3 around the line $\sigma_{1}^{2}=\sigma_{2}^{2}$. Nevertheless, very large noise variances can drown out the impact of small noise variances, regardless of their relative proportions. This behavior provides a rough explanation for the sensitivity of PCA to even a few gross errors (i.e., outliers); even in small proportions, sufficiently large errors dominate the performance of PCA.

Along the dashed cyan line in Figure 3, the average noise variance is $\bar{\sigma}^{2}\approx 1.74$ and the best performance occurs when $\sigma_{1}^{2}=\sigma_{2}^{2}=\bar{\sigma}^{2}$, as predicted by Theorem 2. Along the dotted green curve, the average inverse noise variance is $\mathcal{I}\approx 0.575$ and the best performance again occurs when $\sigma_{1}^{2}=\sigma_{2}^{2}$, as predicted in Remark 6. Note, in particular, that the dashed line and dotted curve are both tangent to the contour at exactly $\sigma_{1}^{2}=\sigma_{2}^{2}$. The observation that larger noise variances have “more impact” provides a rough explanation for this phenomenon; homoscedasticity minimizes the largest noise variance for both the line and the curve. In some sense, as discussed in Section 2.6, the degradation from samples with larger noise is greater than the benefit of having samples with correspondingly smaller noise.

3.4 Impact of adding data

Consider adding data with noise variance $\sigma_{2}^{2}$ and sample-to-dimension ratio $c_{2}$ to an existing dataset that has noise variance $\sigma_{1}^{2}=1$, sample-to-dimension ratio $c_{1}=10$ and subspace amplitude $\theta_{i}=1$ for the $i$th component. The combined dataset has a sample-to-dimension ratio of $c=c_{1}+c_{2}$ and is potentially heteroscedastic with noise variances $\sigma_{1}^{2}$ and $\sigma_{2}^{2}$ appearing in proportions $p_{1}=c_{1}/c$ and $p_{2}=c_{2}/c$. Figure 4 shows the asymptotic subspace recovery (4) of the $i$th component for this combined dataset as a function of the sample-to-dimension ratio $c_{2}$ of the added data for a variety of noise variances $\sigma_{2}^{2}$. The dashed orange curve, showing the recovery when $\sigma_{2}^{2}=1=\sigma_{1}^{2}$, illustrates the benefit we would expect for homoscedastic data: increasing the samples per dimension improves recovery. The dotted red curve shows the recovery when $\sigma_{2}^{2}=4>\sigma_{1}^{2}$. For a small number of added samples, the harm of introducing noisier data outweighs the benefit of having more samples. For sufficiently many samples, however, the tradeoff reverses and recovery for the combined dataset exceeds that for the original dataset; the break even point can be calculated using expression (4). Finally, the green curve shows the performance when $\sigma_{2}^{2}=1.4>\sigma_{1}^{2}$. As before, the added samples are noisier than the original samples and so we might expect performance to initially decline again. In this case, however, the performance improves for any number of added samples. In all three cases, the added samples dominate in the limit $c_{2}\to\infty$ and PCA approaches perfect subspace recovery as one may expect. However, perfect recovery in the limit does not typically happen for PCA amplitudes (2) and coefficient recovery (5); see Section S3.4 of the Online Supplement for more details.

Note that it is equivalent to think about removing noisy samples from a dataset by thinking of the combined dataset as the original full dataset. The green curve in Figure 4 then suggests that slightly noisier samples should not be removed; it would be best if the full data was homoscedastic but removing slightly noisier data (and reducing the dataset size) does more harm than good. The dotted red curve in Figure 4 suggests that much noisier samples should be removed unless they are numerous enough to outweigh the cost of adding them. Once again, expression (4) can be used to calculate the break even point.

4 Numerical simulation

This section simulates data generated by the model described in Section 2.1 to illustrate the main result, Theorem 1, and to demonstrate that the asymptotic results provided are meaningful for practical settings with finitely many samples in a finite-dimensional space. As in Section 3, we show results only for the asymptotic subspace recovery (4) for brevity; the same phenomena occur for the asymptotic PCA amplitudes (2) and coefficient recovery (5) as we show in Section S4 of the Online Supplement. Consider data from a two-dimensional subspace with subspace amplitudes $\theta_{1}=1$ and $\theta_{2}=0.8$, two noise variances $\sigma_{1}^{2}=0.1$ and $\sigma_{2}^{2}=3.25$, and a sample-to-dimension ratio of $c=10$. We sweep the proportion of high noise samples $p_{2}$ from zero to one, setting $p_{1}=1-p_{2}$ as in Section 3.2. The first simulation considers $n=10^{3}$ samples in a $d=10^{2}$ dimensional ambient space ($10^{4}$ trials). The second increases these to $n=10^{4}$ samples in a $d=10^{3}$ dimensional ambient space ($10^{3}$ trials). Both simulations generate data from the standard normal distribution, i.e., $z_{ij},\varepsilon_{ij}\sim\mathcal{N}(0,1)$. Note that sweeping over $p_{2}$ covers homoscedastic settings at the extremes ($p_{2}=0,1$) and evenly split heteroscedastic data in the middle ($p_{2}=1/2$).

Figure 5 plots the recovery of subspace components $|\langle\hat{u}_{i},u_{i}\rangle|^{2}$ for both simulations with the mean (dashed blue curve) and interquartile interval (light blue ribbon) shown with the asymptotic recovery (4) of Theorem 1 (green curve). The region where $A(\beta_{i})\leq 0$ is the red horizontal segment with value zero (the prediction of Conjecture 1). Figure 5(a) illustrates general agreement between the mean and the asymptotic recovery, especially far away from the non-differentiable points where the recovery becomes zero and Conjecture 1 predicts a phase transition. This is a general phenomenon we observed: near the phase transition the smooth simulation mean deviates from the non-smooth asymptotic recovery. Intuitively, an asymptotic recovery of zero corresponds to PCA components that are like isotropically random vectors and so have vanishing square inner product with the true components as the dimension grows. In finite dimension, however, there is a chance of alignment that results in a positive square inner product.

Figure 5(b) shows what happens when the number of samples and ambient dimension are increased to $n=10^{4}$ and $d=10^{3}$. The interquartile intervals are roughly half the size of those in Figure 5(a), indicating concentration of the recovery of each component (a random quantity) around its mean. Furthermore, there is better agreement between the mean and the asymptotic recovery, with the maximum deviation between simulation and asymptotic prediction still occurring nearby the phase transition. In particular for $p_{2}<0.75$ the largest deviation for $|\langle\hat{u}_{1},u_{1}\rangle|^{2}$ is around $0.03$. For $p_{2}\notin(0.1,0.35)$, the largest deviation for $|\langle\hat{u}_{2},u_{2}\rangle|^{2}$ is around $0.02$. To summarize, the numerical simulations indicate that the subspace recovery concentrates to its mean and that the mean approaches the asymptotic recovery. Furthermore, good agreement with Conjecture 1 provides further evidence that there is indeed a phase transition below which the subspace is not recovered. These findings are similar to those in [18] for a one-dimensional subspace with two noise variances.

5 Proof of Theorem 1

The proof has six main parts. Section 5.1 connects several results from random matrix theory to obtain an initial expression for asymptotic recovery. This expression is difficult to evaluate and analyze because it involves an integral transform of the (nontrivial) limiting singular value distribution for a random (noise) matrix as well as the corresponding limiting largest singular value. However, we have discovered a nontrivial structure in this expression that enables us to derive a much simpler form in Sections 5.2-5.6.

5.1 Obtain an initial expression

Rewriting the model in (1) in matrix form yields

[TABLE]

where

• $\mathbf{Z}=(z^{(1)},\ldots,z^{(k)})\in\mathbb{C}^{n\times k}$ is the coefficient matrix,
• $\mathbf{E}=(\varepsilon_{1},\ldots,\varepsilon_{n})\in\mathbb{C}^{d\times n}$ is the (unscaled) noise matrix,
• $\mathbf{H}=\operatorname*{diag}(\eta_{1},\dots,\eta_{n})\in\mathbb{R}_{+}^{n\times n}$ is a diagonal matrix of noise standard deviations.

The first $k$ principal components $\hat{u}_{1},\ldots,\hat{u}_{k}$, PCA amplitudes $\hat{\theta}_{1},\ldots,\hat{\theta}_{k}$ and (normalized) scores $\hat{z}^{(1)}/\sqrt{n},\ldots,\hat{z}^{(k)}/\sqrt{n}$ defined in Section 1 are exactly the first $k$ left singular vectors, singular values and right singular vectors, respectively, of the scaled data matrix $\mathbf{Y}/\sqrt{n}$.

To match the model of [5], we introduce the random unitary matrix

[TABLE]

where the random matrix $\mathbf{\tilde{U}}\in\mathbb{C}^{d\times k}$ is the Gram-Schmidt orthonormalization of a $d\times k$ random matrix that has iid (mean zero, variance one) circularly symmetric complex normal $\mathcal{CN}(0,1)$ entries. We use the superscript $\perp$ to denote a matrix of orthonormal basis elements for the orthogonal complement; the columns of $\mathbf{U}^{\perp}$ form an orthonormal basis for the orthogonal complement of the column span of $\mathbf{U}$.

Left multiplying (17) by $\mathbf{R}/\sqrt{n}$ yields that $\mathbf{R}\hat{u}_{1},\ldots,\mathbf{R}\hat{u}_{k}$, $\hat{\theta}_{1},\ldots,\hat{\theta}_{k}$ and $\hat{z}^{(1)}/\sqrt{n},\ldots,\hat{z}^{(k)}/\sqrt{n}$ are the first $k$ left singular vectors, singular values and right singular vectors, respectively, of the scaled and rotated data matrix

[TABLE]

The matrix $\mathbf{\tilde{Y}}$ matches the low rank (i.e., rank $k$) perturbation of a random matrix model considered in [5] because

[TABLE]

where

[TABLE]

Here $\mathbf{P}$ is generated according to the “orthonormalized model” in [5] for the vectors $\tilde{u}_{i}$ and the “iid model” for the vectors $z^{\left(i\right)}$ and $\mathbf{P}$ satisfies Assumption 2.4 of [5]; the latter considers $\tilde{u}_{i}$ and $z^{\left(i\right)}$ to be generated according to the same model, but its proof extends to this case. Furthermore $\mathbf{R}\mathbf{E}$ has iid entries with zero mean, unit variance and bounded fourth moment (by the assumption that $\varepsilon_{i}$ are unitarily invariant), and $\mathbf{H}$ is a non-random diagonal positive definite matrix with bounded spectral norm and limiting eigenvalue distribution $p_{1}\delta_{\sigma_{1}^{2}}+\cdots+p_{L}\delta_{\sigma_{L}^{2}}$, where $\delta_{\sigma_{\ell}^{2}}$ is the Dirac delta distribution centered at $\sigma_{\ell}^{2}$. Under these conditions, Theorem 4.3 and Corollary 6.6 of [3] state that $\mathbf{X}$ has a non-random compactly supported limiting singular value distribution $\mu_{\mathbf{X}}$ and the largest singular value of $\mathbf{X}$ converges almost surely to the supremum of the support of $\mu_{\mathbf{X}}$. Thus Assumptions 2.1 and 2.3 of [5] are also satisfied.

Furthermore, $\hat{u}_{i}^{\mathsf{H}}u_{j}=\hat{u}_{i}^{\mathsf{H}}\mathbf{R}^{\mathsf{H}}\mathbf{R}u_{j}=(\mathbf{R}\hat{u}_{i})^{\mathsf{H}}\tilde{u}_{j}$ for all $i,j\in\{1,\dots,k\}$ so

[TABLE]

and hence Theorem 2.10 from [5] implies that, for each $i\in\{1,\ldots,k\}$,

[TABLE]

and that if $\theta_{i}^{2}>\bar{\theta}^{2}$, then

[TABLE]

and

[TABLE]

where $\rho_{i}=D^{-1}(1/\theta_{i}^{2})$, $\bar{\theta}^{2}=1/D\left(b^{+}\right)$, $D(z)=\varphi(z)\{c^{-1}\varphi\left(z\right)+(1-c^{-1})/z\}$ for $z>b$, $\varphi\left(z\right)=\int z/(z^{2}-t^{2})\ d\mu_{\mathbf{X}}\left(t\right),$ $b$ is the supremum of the support of $\mu_{\mathbf{X}}$ and $\mu_{\mathbf{X}}$ is the limiting singular value distribution of $\mathbf{X}$ (compactly supported by Assumption 2.1 of [5]). We use the notation $f\left(b^{+}\right)=\lim_{z\rightarrow b^{+}}f\left(z\right)$ as a convenient shorthand for the limit from above of a function $f\left(z\right)$.

Theorem 2.10 from [5] is presented therein for $d\leq n$ (i.e., $c\geq 1$) to simplify their proofs. However, it also holds without modification for $d>n$ if the limiting singular value distribution $\mu_{\mathbf{X}}$ is always taken to be the limit of the empirical distribution of the $d$ largest singular values ($d-n$ of which will be zero). Thus we proceed without the condition that $c>1$.

Furthermore, even though it is not explicitly stated as a main result in [5], the proof of Theorem 2.10 in [5] implies that

[TABLE]

as was also noted in [27] for the special case of distinct subspace amplitudes.

Evaluating the expressions (18), (19) and (21) would consist of evaluating the intermediates listed above from last to first. These steps are challenging because they involve an integral transform of the limiting singular value distribution $\mu_{\mathbf{X}}$ for the random (noise) matrix $\mathbf{X}$ as well as the corresponding limiting largest singular value $b$, both of which depend nontrivially on the model parameters. Our analysis uncovers a nontrivial structure that we exploit to derive simpler expressions.

Before proceeding, observe that the almost sure limit in (21) is just the geometric mean of the two almost sure limits in (19). Hence, we proceed to derive simplified expressions for (18) and (19); (6) follows as the geometric mean of the simplified expressions obtained for the almost sure limits in (19).

5.2 Perform a change of variables

We introduce the function defined, for $z>b$, by

[TABLE]

because it turns out to have several nice properties that simplify all of the following analysis. Rewriting (19) using $\psi\left(z\right)$ instead of $\varphi\left(z\right)$ and factoring appropriately yields that if $\theta_{i}^{2}>\bar{\theta}^{2}$ then

[TABLE]

where now

[TABLE]

for $z>b$ and we have used the fact that

[TABLE]

5.3 Find useful properties of $\psi(z)$

Establishing some properties of $\psi\left(z\right)$ aids simplification significantly.

Property 1. We show that $\psi\left(z\right)$ satisfies a certain rational equation for all $z>b$ and derive its inverse function $\psi^{-1}(x)$. Observe that the square singular values of the noise matrix $\mathbf{X}$ are exactly the eigenvalues of $c\mathbf{X}\mathbf{X}^{\mathsf{H}}$, divided by $c$. Thus we first consider the limiting eigenvalue distribution $\mu_{c\mathbf{X}\mathbf{X}^{\mathsf{H}}}$ of $c\mathbf{X}\mathbf{X}^{\mathsf{H}}$ and then relate its Stieltjes transform $m\left(\zeta\right)$ to $\psi\left(z\right)$.

Theorem 4.3 in [3] establishes that the random matrix $c\mathbf{X}\mathbf{X}^{\mathsf{H}}=(1/d)\mathbf{E}\mathbf{H}^{2}\mathbf{E}^{\mathsf{H}}$ has a limiting eigenvalue distribution $\mu_{c\mathbf{X}\mathbf{X}^{\mathsf{H}}}$ whose Stieltjes transform is given, for $\zeta\in\mathbb{C}^{+}$, by

[TABLE]

and satisfies the condition

[TABLE]

where $\mathbb{C}^{+}$ is the set of all complex numbers with positive imaginary part.

Since the $d$ square singular values of $\mathbf{X}$ are exactly the $d$ eigenvalues of $c\mathbf{X}\mathbf{X}^{\mathsf{H}}$ divided by $c$, we have for all $z>b$

[TABLE]

For all $z$ and $\xi>0$, $z^{2}c+i\xi\in\mathbb{C}^{+}$ and so combining (25)–(27) yields that for all $z>b$

[TABLE]

Rearranging yields

[TABLE]

for all $z>b$, where the last term is

[TABLE]

because $p_{1}+\cdots+p_{L}=1$. Substituting back into (28) yields $0=Q\{\psi\left(z\right),z\}$ for all $z>b$, where

[TABLE]

Thus $\psi(z)$ is an algebraic function (the associated polynomial can be formed by clearing the denominator of $Q$). Solving (29) for $z>b$ yields the inverse

[TABLE]

Property 2. We show that $\max_{\ell}(\sigma_{\ell}^{2})<\psi(z)<cz^{2}$ for $z>b$. For $z>b$, one can show from (22) that $\psi(y)$ increases continuously and monotonically from $\psi(z)$ to infinity as $y$ increases from $z$ to infinity, and hence $\psi^{-1}(x)$ must increase continuously and monotonically from $z$ to infinity as $x$ increases from $\psi(z)$ to infinity. However, $\psi^{-1}(x)$ is discontinuous at $x=\max_{\ell}(\sigma_{\ell}^{2})$ because $\psi^{-1}(x)\to\infty$ as $x\to\max_{\ell}(\sigma_{\ell}^{2})$ from the right, and so it follows that $\psi(z)>\max_{\ell}(\sigma_{\ell}^{2})$. Thus $1/\{\psi(z)-\sigma_{\ell}^{2}\}>0$ for all $\ell\in\{1,\dots,L\}$ and so

[TABLE]

Property 3. We show that $0<\psi\left(b^{+}\right)<\infty$ and $\psi^{\prime}\left(b^{+}\right)=\infty$. Property 2 in the limit $z=b^{+}$ implies that

[TABLE]

Taking the total derivative of $0=Q\{\psi(z),z\}$ with respect to $z$ and solving for $\psi^{\prime}(z)$ yields

[TABLE]

As observed in [28], the non-pole boundary points of compactly supported distributions like $\mu_{c\mathbf{X}\mathbf{X}^{\mathsf{H}}}$ occur where the polynomial defining their Stieltjes transform has multiple roots. Thus $\psi(b^{+})$ is a multiple root of $Q(\cdot,b)$ and so

[TABLE]

Thus $\psi^{\prime}\left(b^{+}\right)=\infty$, where the sign is positive because $\psi\left(z\right)$ is an increasing function on $z>b$.

Summarizing, we have shown that

1. a)

$0=Q\{\psi\left(z\right),z\}$ for all $z>b$ where $Q$ is defined in (29), and the inverse function $\psi^{-1}(x)$ is given in (30), 2. b)

$\max_{\ell}(\sigma_{\ell}^{2})<\psi(z)<cz^{2}$, 3. c)

$0<\psi\left(b^{+}\right)<\infty$ and $\psi^{\prime}\left(b^{+}\right)=\infty$.

We now use these properties to aid simplification.

5.4 Express $D(z)$ and $D^{\prime}(z)/z$ in terms of only $\psi(z)$

We can rewrite (24) as

[TABLE]

because $0=Q\{\psi\left(z\right),z\}$ by Property 1 of Section 5.3. Differentiating (32) with respect to $z$ yields

[TABLE]

and so we need to find $\psi^{\prime}\left(z\right)$ in terms of $\psi\left(z\right)$. Substituting the expressions for the partial derivatives $\partial Q\{\psi\left(z\right),z\}/\partial z$ and $\partial Q\{\psi\left(z\right),z\}/\partial s$ into (31) and simplifying we obtain $\psi^{\prime}\left(z\right)=2cz/\gamma\left(z\right)$, where the denominator is

[TABLE]

Note that

[TABLE]

because $0=Q\{\psi\left(z\right),z\}$ for $z>b$. Substituting into $\gamma(z)$ and forming a common denominator, then dividing with respect to $\psi(z)$ yields

[TABLE]

where $A(x)$ was defined in (3). Thus

[TABLE]

and

[TABLE]

where $B_{i}^{\prime}(x)$ is the derivative of $B_{i}(x)$ defined in (3).

5.5 Express the asymptotic recoveries in terms of only $\psi(b^{+})$ and $\psi(\rho_{i})$

Evaluating (32) in the limit $z=b^{+}$ and recalling that $D(b^{+})=1/\bar{\theta}^{2}$ yields

[TABLE]

where $B_{i}(x)$ was defined in (3). Evaluating the inverse function (30) both for $\psi(\rho_{i})$ and in the limit $\psi(b^{+})$ then substituting into (18) yields

[TABLE]

Evaluating (34) for $z=\rho_{i}$ and substituting into (23) yields

[TABLE]

if $B_{i}\{\psi\left(b^{+}\right)\}<0$.

5.6 Obtain algebraic descriptions

This subsection obtains algebraic descriptions of (35), (36) and (37) by showing that $\psi(b^{+})$ is the largest real root of $A(x)$ and that $\psi(\rho_{i})$ is the largest real root of $B_{i}(x)$ when $\theta_{i}^{2}>\bar{\theta}^{2}$. Evaluating (33) in the limit $z=b^{+}$ yields

[TABLE]

because $\psi^{\prime}\left(b^{+}\right)=\infty$ by Property 3 of Section 5.3. If $\theta_{i}^{2}>\bar{\theta}^{2}$ then $\rho_{i}=D^{-1}(1/\theta_{i}^{2})$ and so

[TABLE]

(38) shows that $\psi\left(b^{+}\right)$ is a real root of $A(x)$, and (39) shows that $\psi\left(\rho_{i}\right)$ is a real root of $B_{i}(x)$.

Recall that $\psi(b^{+}),\psi(\rho_{i})\geq\max_{\ell}(\sigma_{\ell}^{2})$ by Property 2 of Section 5.3, and note that both $A(x)$ and $B_{i}(x)$ monotonically increase for $x>\max_{\ell}(\sigma_{\ell}^{2})$. Thus each has exactly one real root larger than $\max_{\ell}(\sigma_{\ell}^{2})$, i.e., its largest real root, and so $\psi(b^{+})=\alpha$ and $\psi(\rho_{i})=\beta_{i}$ when $\theta_{i}^{2}>\bar{\theta}^{2}$, where $\alpha$ and $\beta_{i}$ are the largest real roots of $A(x)$ and $B_{i}(x)$, respectively.

A subtle point is that $A(x)$ and $B_{i}(x)$ always have largest real roots $\alpha$ and $\beta$ even though $\psi(\rho_{i})$ is defined only when $\theta_{i}^{2}>\bar{\theta}^{2}$. Furthermore, $\alpha$ and $\beta$ are always larger than $\max_{\ell}(\sigma_{\ell}^{2})$ and both $A(x)$ and $B_{i}(x)$ are monotonically increasing in this regime and so we have the equivalence

[TABLE]

Writing (35), (36) and (37) in terms of $\alpha$ and $\beta_{i}$, then applying the equivalence (40) and combining with (20) yields the main results (2), (4) and (5).

6 Proof of Theorem 2

If $A(\beta_{i})\geq 0$ then (4) and (5) increase with $A(\beta_{i})$ and decrease with $\beta_{i}$ and $B^{\prime}(\beta_{i})$. Similarly, (2) increases with $\beta_{i}$, as illustrated by (5). As a result, Theorem 2 follows immediately from the following bounds, all of which are met with equality if and only if $\sigma_{1}^{2}=\cdots=\sigma_{L}^{2}$:

[TABLE]

The bounds (41) are shown by exploiting convexity to appropriately bound the rational functions $B_{i}(x)$, $B_{i}^{\prime}(x)$ and $A(x)$. We bound $\beta_{i}$ by noting that

[TABLE]

because $\sigma_{\ell}^{2}<\beta_{i}$ and $f(v)=1/(\beta_{i}-v)$ is a strictly convex function over $v<\beta_{i}$. Thus $\beta_{i}\geq c\theta_{i}^{2}+\bar{\sigma}^{2}$. We bound $B_{i}^{\prime}(\beta_{i})$ by noting that

[TABLE]

because the quadratic function $z^{2}$ is strictly convex. Similarly,

[TABLE]

because the quadratic function $z^{2}$ is strictly convex and

[TABLE]

All of the above bounds are met with equality if and only if $\sigma_{1}^{2}=\cdots=\sigma_{L}^{2}$ because the convexity in all cases is strict. As a result, homoscedastic noise minimizes (2), and it maximizes (4) and (5). See Section S2 of the Online Supplement for some interesting additional properties in this context.

7 Discussion and extensions

This paper provided simplified expressions (Theorem 1) for the asymptotic recovery of a low-dimensional subspace, the corresponding subspace amplitudes and the corresponding coefficients by the principal components, PCA amplitudes and scores, respectively, obtained from applying PCA to noisy high-dimensional heteroscedastic data. The simplified expressions provide generalizations of previous results for the special case of homoscedastic data. They were derived by first connecting several recent results from random matrix theory [3, 5] to obtain initial expressions for asymptotic recovery that are difficult to evaluate and analyze, and then exploiting a nontrivial structure in the expressions to find the much simpler algebraic descriptions of Theorem 1.

These descriptions enable both easy and efficient calculation as well as reasoning about the asymptotic performance of PCA. In particular, we use the simplified expressions to show that, for a fixed average noise variance, asymptotic subspace recovery, amplitude recovery and coefficient recovery are all worse when the noise is heteroscedastic as opposed to homoscedastic (Theorem 2). Hence, while average noise variance is often a practically convenient measure for the overall quality of data, it gives an overly optimistic estimate of PCA performance. Our expressions (2), (4) and (5) in Theorem 1 are more accurate.

We also investigated examples to gain insight into how the asymptotic performance of PCA depends on the model parameters: sample-to-dimension ratio $c$, subspace amplitudes $\theta_{1},\dots,\theta_{k}$, proportions $p_{1},\dots,p_{L}$ and noise variances $\sigma_{1}^{2},\dots,\sigma_{L}^{2}$. We found that performance depends in expected ways on

• a)

sample-to-dimension ratio: performance improves with more samples;

• b)

subspace amplitudes: performance improves with larger amplitudes;

• c)

proportions: performance improves when more samples have low noise.

We also learned that when the gap between the two largest noise variances is “sufficiently wide”, the performance is dominated by the largest noise variance. This result provides insight into why PCA performs poorly in the presence of gross errors and why heteroscedasticity degrades performance in the sense of Theorem 2. Nevertheless, adding “slightly” noisier samples to an existing dataset can still improve PCA performance; even adding significantly noisier samples can be beneficial if they are sufficiently numerous.

Finally, we presented numerical simulations that demonstrated concentration of subspace recovery to the asymptotic prediction (4) with good agreement for practical problem sizes. The same agreement occurs for the PCA amplitudes and coefficient recovery. The simulations also showed good agreement with the conjectured phase transition (Conjecture 1).

There are many exciting avenues for extensions and further work. An area of ongoing work is the extension of our analysis to a weighted version of PCA, where the samples are first weighted to reduce the impact of noisier points. Such a method may be natural when the noise variances are known or can be estimated well. Data formed in this way do not match the model of [5], and so the analysis involves first extending the results of [5] to handle this more general case. Preliminary findings suggest that whitening the noise with inverse noise variance weights $1/\sigma_{\ell}^{2}$ is not optimal.

Another natural direction is to consider a general distribution of noise variances $\nu$, where we suppose that the empirical noise distribution $(\delta_{\eta_{1}^{2}}+\cdots+\delta_{\eta_{n}^{2}})/n\overset{a.s.}{\longrightarrow}\nu$ as $n\to\infty$. We conjecture that if $\eta_{1},\dots,\eta_{n}$ are bounded for all $n$ and $\int d\nu(\tau)/(x-\tau)\to\infty$ as $x\to\tau_{\text{max}}^{+}$, then the almost sure limits in this paper hold but with

[TABLE]

where $\tau_{\text{max}}$ is the supremum of the support of $\nu$. The proofs of Theorem 1 and Theorem 2 both generalize straightforwardly for the most part; the main trickiness comes in carefully arguing that limits pass through integrals in Section 5.3.

Proving that there is indeed a phase transition in the asymptotic subspace recovery and coefficient recovery, as conjectured in Conjecture 1, is another area of future work. That proof may be of greater interest in the context of a weighted PCA method. Another area of future work is explaining the puzzling phenomenon described in Section 3.3, where, in some regimes, performance improves by increasing the noise variance. More detailed analysis of the general impacts of the model parameters could also be interesting. A final direction of future work is deriving finite sample results for heteroscedastic noise as was done for homoscedastic noise in [29].

Acknowledgments

The authors thank Raj Rao Nadakuditi and Raj Tejas Suryaprakash for many helpful discussions regarding the singular values and vectors of low rank perturbations of large random matrices. The authors also thank Edgar Dobriban for his feedback on a draft and for pointing them to the generalized spiked covariance model. The authors also thank Rina Foygel Barber for suggesting average inverse noise variance be considered and for pointing out that Theorem 2 implies an analogous statement for this measure. Finally, the authors thank the editors and referees for their many helpful comments and suggestions that significantly improved the strength, clarity and style of the paper.