Resampling Sensitivity of High-Dimensional PCA

Let $\mathbf{X}=(\mathbf{X}_{ij})$ be an $n\times p$ data matrix with independent real valued entries with mean 0 and variance $p^{-1}$,

Note that we do not require the i.i.d. condition for the data. Furthermore, we assume the entries $x_{ij}$ have a sub-exponential decay, that is, there exists a constant $\theta>0$ such that for $u>1$,

This sub-exponential decay assumption is mainly for convenience, and other conditions such as the finiteness of a sufficiently high moment would be enough.

Motivated by high-dimensional statistics, we will work in the proportional growth regime $n\asymp p$. Throughout this paper, to avoid trivial eigenvalues, we will be working in the regime

In the case $\lim p/n=1$, our assymption $p/n\equiv 1$ is due to technical reasons in random matrix theory. Specifically, the proof relies on the delocalization of eigenvectors in the whole spectrum. As one of the major open problems in random matrix theory, delocalization of eigenvectors near the lower spectral edge is not known in the general case with just $\lim p/n=1$. The strictly square assumption $p\equiv n$ can be slightly relaxed to $|n-p|=p^{o(1)}$ (see e.g. [Wan22]), but we do not pursue such an extension for simplicity.

The sample covariance matrix corresponding to data matrix $\mathbf{X}$ is defined by $\mathbf{H}:=\mathbf{X}^{\top}\mathbf{X}$. We order the eigenvalues of $\mathbf{H}$ as $\lambda_{1}\geq\cdots\geq\lambda_{p}$, and use $\mathbf{v}_{i}\in\mathbb{R}^{p}$ to denote the unit eigenvector corresponding to the eigenvalue $\lambda_{i}$. If the context is clear, we just use $\lambda:=\lambda_{1}$ and $\mathbf{v}:=\mathbf{v}_{1}$ to denote the largest eigenvalue and the top eigenvector. We also consider the eigenvalues and eigenvectors of the Gram matrix $\widehat{\mathbf{H}}:=\mathbf{X}\mathbf{X}^{\top}$. Note that $\widehat{\mathbf{H}}$ and $\mathbf{H}$ have the same non-trivial eigenvalues, and the spectrum of $\widehat{\mathbf{H}}$ is given by $\{\lambda_{i}\}_{i=1}^{n}$ with $\lambda_{p+1}=\cdots=\lambda_{n}=0$. We denote the unit eigenvectors of $\widehat{\mathbf{H}}$ associated with the eigenvalue $\lambda_{i}$ by $\mathbf{u}_{i}\in\mathbb{R}^{n}$.

Let $\mathbf{U}=[\mathbf{u}_{1},\cdots,\mathbf{u}_{p}]\in\mathbb{R}^{n\times p}$ and $\mathbf{V}=[\mathbf{v}_{1},\cdots,\mathbf{v}_{p}]\in\mathbb{R}^{p\times p}$. These eigenvectors may be connected by the singular value decomposition of the data matrix $\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$, where $\mathbf{\Sigma}:=\mathsf{diag}(\sigma_{1},\cdots,\sigma_{p})$ with $\sigma_{i}=\sqrt{\lambda_{i}}$ corresponds to the singular values. For convenience, we also denote $\sigma:=\sigma_{1}$. And therefore, up to the sign of the eigenvectors, we have

We now describe the resampling procedure. For a positive number $k\leq np$, define the random matrix $\mathbf{X}^{[k]}$ in the following way. Let $S_{k}=\{(i_{1},\alpha_{1}),\cdots,(i_{k},\alpha_{k})\}$ be a set of $k$ pairs chosen uniformly at random without raplacement from the set of all ordered pairs $(i,\alpha)$ of indices with $1\leq i\leq n$ and $1\leq\alpha\leq p$. We assume that the set $S_{k}$ is independent of the entries of $\mathbf{X}$. The entries of $\mathbf{X}^{[k]}$ are given by

where $(\mathbf{X}_{i,\alpha}^{\prime})_{1\leq i\leq n,1\leq\alpha\leq p}$ are independent random variables, independent of $\mathbf{X}$, and $\mathbf{X}_{i,\alpha}^{\prime}$ has the same distribution as $\mathbf{X}_{i,\alpha}$. In other words, $\mathbf{X}^{[k]}$ is obtained from $\mathbf{X}$ by resampling $k$ random entries of the matrix, and therefore $\mathbf{X}^{[k]}$ clearly has the same distribution as $\mathbf{X}$. Let $\mathbf{H}^{[k]}:=(\mathbf{X}^{[k]})^{\top}\mathbf{X}^{[k]}$ the sample covariance matrix corresponding to the resampled data matrix $\mathbf{X}^{[k]}$. Denote the eigenvalues and the corresponding normalized eigenvectors of $\mathbf{H}^{[k]}$ by $\lambda^{[k]}_{1}\geq\cdots\geq\lambda^{[k]}_{p}$ and $\mathbf{v}^{[k]}_{1},\cdots,\mathbf{v}^{[k]}_{p}$. When the context is clear, the principal component is just denoted by $\lambda^{[k]}$ and $\mathbf{v}^{[k]}$. Similarly, denote the eigenvector of the matrix $\widehat{\mathbf{H}}^{[k]}:=\mathbf{X}^{[k]}(\mathbf{X}^{[k]})^{\top}$ associated with the eigenvalue $\lambda^{[k]}_{i}$ by $\mathbf{u}^{[k]}_{i}$.

With the resampling parameter in two different regimes, we have the following results.

These two theorems together state that the critical threshold for the resampling strength is of order $k\asymp n^{5/3}$. Note that $n^{5/3}$ compared with the total number of inputs $np\asymp n^{2}$ is negligible. We show that, above the threshold $n^{5/3}$, resampling even a negligible portion of the data will result in a dramatic change of the resulting principal component, in the sense that the new principal component is asymptotically orthogonal to the old one; while below the threshold, a relatively mild resampling has almost no effect on the corresponding new principal component. If considering the eigenvector overlaps $|\langle\mathbf{v},\mathbf{v}^{[k]}\rangle|$ and $|\langle\mathbf{u},\mathbf{u}^{[k]}\rangle|$, these quantities exhibit sharp phase transitions from $1$ to $0$ near the critical threshold $k\asymp n^{5/3}$.

We remark that the phase transition stated in the above theorems is not restricted to the top eigenvectors $\mathbf{v},\mathbf{v}^{[k]},\mathbf{u},\mathbf{u}^{[k]}$. With essentially the same arguments, we can prove that for any fixed $m>0$, the $m$-th leading eigenvectors $\mathbf{v}_{m},\mathbf{v}^{[k]}_{m}$ and $\mathbf{u}_{m},\mathbf{u}^{[k]}_{m}$ exhibit the same phase transition at the critical threshold of the same order $n^{5/3}$.

The high-level idea is the “superconcentration implies chaos” phenomenon established by Chatterjee [Cha14], which means that small perturbation (beyond certain threshold) of super-concentrated system ( in the sense that it is characterized by some quantity with small variance) will lead to a dramatic change of the system. For random matrix models, the super-concentrated quantity is usually the eigenvalues. Resampling sensitivity for random matrices was first studied for Wigner matrices and sparse Erdős-Rényi graphs in [BLZ20, BL22]. For the PCA model, the key difference is that the entries of the sample covariance matrix are correlated. Moreover, resampling a single entry of the data will change $\Theta(n)$ entries in the sample covariance matrix. These two differences will make the proofs more technical and a linearization trick would be important to reduce the interdependency of the matrix entries. In addition, we also need a technical variance formula related to resampling to compute the variance of the top eigenvalue and tools from random matrix theory, in particular the local Marchenko-Pastur law for the resolvent and the delocalization of eigenvectors.

For the sensitivity of the principal components. We show that the inner products $\langle\mathbf{v},\mathbf{v}^{[k]}\rangle$ and $\langle\mathbf{u},\mathbf{u}^{[k]}\rangle$ are closely related to the variance of the top eigenvalue of the sample covariance matrix. Using the variance formula for resampling, we have a precise characterization between the inner product of the top eigenvectors and the concentration of the top eigenvalue. The perturbation effect of the resampling is studied via the variational representation of the eigenvalue.

For the stability under moderate resampling, the key idea of our proof is to study the stability of the resolvent. On the one hand, the stability of resolvent implies the stability of the top eigenvalue. On the other hand, the resolvent can be used to approximate certain useful statistics of the top eigenvector. The stability of the resolvent is proved via a Lindeberg exchange strategy. The resampling procedure can be decomposed into a martingale, and the difference between the resolvents can be therefore bounded by martingale concentration inequalities combined with the local Marchenko-Pastur law of the resolvent for a priori estimates.

As mentioned above, the key step in the proof for sensitivity regime is to establish a connection between the inner products of top eigenvalues with the variance of the top eigenvalue. Specifically, we will prove

where $C>0$ is some universal constant and $\sigma=\sqrt{\lambda}$ is the leading singular value. A similar result is also true for $\mathbf{u}$ and $\mathbf{u}^{[k]}$. For more details, we refer to Section C.2.

From random matrix theory, we have $\mathsf{Var}(\sigma)=O(n^{-4/3})$. Then, based on this inequality, we derive the threshold $k\gg n^{5/3}$ for the sensitivity regime.

To illustrate the usefulness of the resolvent, we show that the entries of the resolvent can be used to approximate the product of entries in the eigenvector. For some small $\delta>0$, let $z_{0}=\lambda+\mathrm{i}\eta$ with $\eta=n^{-2/3-\delta}$. In the regime $k\leq n^{5/3-\epsilon_{0}}$ for some $\epsilon_{0}>0$, there exists some small $c>0$ such that for all $\alpha,\beta=1,\dots,p$, we have

and

A similar result also holds for $\mathbf{u}$ and $\mathbf{u}^{[k]}$. For more details, we refer to Lemma D.5.

Since the eigenvector can be approximated by the resolvent, it suffices to show the stability of the resolvent. Consider the regime $k\leq n^{5/3-\epsilon_{0}}$ for some $\epsilon_{0}>0$. For some small $\delta>0$ and all $z=E+\mathrm{i}\eta$ that is close to the upper spectral edge and $\eta=n^{-2/3-\delta}$, there exists a small constant $c>0$ such that the following is true for all $i,j=1,\dots,n$ and $\alpha,\beta=1,\dots,p$,

and

This is the main technical part of the whole argument, and its proof relies on the Lindeberg exchange method and a martingale concentration argument. For more details, we refer to Lemma D.3.

Combining the stability of the resolvents with the resolvent approximation for eigenvectors, we can conclude that $\mathbf{v}$ and $\mathbf{v}^{[k]}$ must be close (similarly, also for $\mathbf{u}$ and $\mathbf{u}^{[k]}$).

As a byproduct, we derive an stability estimate of the top eigenvalues, which may be of independent interest. For $k\leq n^{5/3-\epsilon_{0}}$ with some arbitrary $\epsilon_{0}>0$ and an arbitrary $\varepsilon>0$,

where $\delta>0$ is some small constant. Note that the fluctuation of the top eigenvalue around its typical location is of order $O(n^{-2/3+\varepsilon})$, this result shows that the top eigenvalues under moderate resampling are indeed non-trivially stable. For more details, we refer to Lemma D.4.

Within the general PCA framework, one important variant is the kernel PCA, which is closely related to the widely used spectral clustering [NJW01, VL07, CBG16]. The corresponding kernel random matrices were studied in [EK10, CS13]. However, the study of these kernel random matrices are far from being well-understood. In particular, the study of eigenvectors were very limited and the indispensable property of delocalization is still unknown.

It would be interesting to explore whether other statistical method share the same resampling sensitivity phenomenon. Random matrices associated with canonical correlation analysis (CCA) or multivariate analysis of variance (MANOVA) are well studied [HPZ16, HPY18, Yan22]. In particular, the Tracy-Widom concentration of the top eigenvalue, one of the most important ingredients that we need, has been proved. We anticipate that these models exhibit a similar stability-sensitivity transition as in PCA.

In our paper, we study the stability of the top eigenvalue and the top eigenvector under resampling in terms of bounding the $\ell_{\infty}$ distance. Such stability estimates can be regarded as the global sensitivity of PCA performed on neighboring datasets. This measurement is closely related to the analysis of differential privacy [DR14]. PCA under differential privacy was previous studied in [BDMN05, CSS13], etc. Our result revisit the problem of designing a private algorithm for solving the principal component. Here we remark that though the statements in Theorem 1 and Theorem 2 are qualitative, a careful examinination of the proof can yield some quantitative estimates. Based on the stability estimates in terms of the $\ell_{\infty}$ metric, a simple Laplace mechanism produces a differentially private version of PCA for computing the top eigenvalue or the top eigenvector. However, compared with [CSS13], our results are limited in the sense that their results are non-asymptotic for all sample size $n$ and data dimension $p$, while ours are restricted to the proportional growth regime.

Moreover, previous works on differentially private PCA focused on neighbouring datasets that differ by one sample vector. Our result may be seen as a refined notion of privacy, since we can analyze the sensitivity of PCA over two “neighbouring” datasets with $k$ different entries for any $k$.

Meanwhile, the largest eigenvalue of the sample covariance matrix plays an important rule in hypothesis testing. For example, the Roy’s largest root test is used in many problems (see e.g. [JN17]). Our result may provide useful insights to construct a differentially private test statistic based on the top eigenvalue.

Database alignment (or in some cases graph matching) refers to the optimization problem in which we are given two datasets and the goal is to find an optimal correspondence between the samples and features that maximally align the data. For datasets $\mathbf{X},\mathbf{Y}\in\mathbb{R}^{n\times p}$, we look for permutations $\pi_{\mathrm{s}}\in{\mathcal{S}}_{n}$ and $\pi_{\mathrm{f}}\in{\mathcal{S}}_{p}$ to solve the optimization problem

where ${\mathcal{S}}_{n}$ and ${\mathcal{S}}_{p}$ are the sets of all permutations on $[n]$ and $[p]$, respectively. This problem is closely related to the Quadratic Assignment Problem (QAP), which is known to be NP-hard to solve or even approximate.

The study of the alignment problem for correlated random databases has a long history. The previous work mainly considered matrices that are correlated through some additive perturbation, and some of the general model were studied with a homogeneous correlation (i.e. the correlation between all correponding pairs are the same). See for example [DCK19, WXS22] and many other works.

Our resampling procedure may be regarded as an adversarial corruption of the dataset, which is a different kind of correlation compared with previous work. To our knowledge, this is the first time to consider database alignment with adversarial corruption. To state the setting of the problem, we have two matrices

where $\mathbf{X}$ is a random matrix satisfying (1) and (2), and $\mathbf{\Pi}_{\mathrm{s}}$ and $\mathbf{\Pi}_{\mathrm{f}}$ are permutation matrices of order $n$ and $p$ chosen uniformly at random. The goal is to recover the permutations $\mathbf{\Pi}_{\mathrm{s}}$ and $\mathbf{\Pi}_{\mathrm{f}}$ based on the observations $\mathbf{X}$ and $\mathbf{Y}$. Here, we can think of $\mathbf{Y}$ as the unlabeled version of $\mathbf{X}$ with adversarial corruption. By considering the covariance matrices, we have

and similarly we consider

A natural idea to reconstruct the permutations $\mathbf{\Pi}_{\mathrm{s}}$ (and $\mathbf{\Pi}_{\mathrm{f}})$ is to align the top eigenvectors of the matrices $\mathbf{A}$ and $\mathbf{B}$ (and $\widehat{\mathbf{A}}$ and $\widehat{\mathbf{B}}$). See Algorithm 1 for details. This spectral method is a natural technique for databse alignment and graph matching. We are interested in under what resampling strength, the PCA-Recovery algorithm can almost perfectly reconstruct the permutations, and under what condition this method completely fail.

Spectral methods have been studied and applied in many scenarios (see e.g. [FMWX20, FMWX22] and [GLM22]). In particular, in [GLM22], a similar PCA method was studied to match two symmetric Gaussian matrices correlated via additive Gaussian noise. Their work proved a similar $0-1$ transition for the inner product of the top eigenvectors, which leads to a all-or-nothing phenomenon in the alignment problem, i.e. the accuracy of the recovery undergoes a sharp transition from $0$ to $1$ near some critical threshold. However, the arguments in [GLM22] are not applicable in our case. Their proof heavily depends on the Gaussian assumption of the matrices, and the additive strucutre of the noise. In particular, they proof crucially relies on the orthogonal invariance of the Gaussian noise. While in our case, the noise is presented in terms of the resampling strength. There is no way to write the “noise” in an additive form that is independent of the “signal”. Even in the Gaussian case, a rigorous analysis of the PCA-Recovery algorithm seems difficult.

Nevertheless, our results on the sensitivity of the eigenvector inner products suggest that, when $k\gg n^{5/3}$, the two eigenvectors are approximately de-correlated so that they share almost no common information. Consequently, recovery of the data via aligning the principal components would be basically random guessing. Therefore, we conjecture that if $k\gg n^{5/3}$, PCA-Recovery fails to recover the latent permutations in the sense that it can only achieve $o(1)$ fraction of correct matching with the ground truth. On the other hand, when $k\ll n^{5/3}$, the performance of our algorithm seems mysterious.

In Section 5, we empirically check the performance of the PCA-Recovery algorithm. Numerical simulations suggest that when $k\gg n^{5/3}$, the performance of PCA-Recovery is indeed poor in the sense that the accuracy of the recovery is almost $0$. On the other hand, when $k\ll n^{5/3}$, experiments show that we cannot expect the sharp all-or-nothing phenomenon similarly as in [GLM22].

Finally, we remark that what PCA-Recovery actually studies is a more difficult task, as we do not need direct observations of $\mathbf{X}$ and $\mathbf{Y}$. We can consider a harder problem (in both statistical and computational sense), which we call alignment from covariance profile. In this problem, we only have access to the covariance between the samples and we aim to recover the correspondence between the samples from the two databases. A similar problem with Gaussian data and additive noise was considered in [WWXY22] as a prototype for matching random geometric dot-product graphs. The analysis of such database alignment problem with adversarial corruption will be an interesting direction for future studies.