Phase Transition in the Recovery of Rank One Matrices Corrupted by Gaussian Noise

TL;DR

This paper investigates the phase transition phenomenon in recovering rank-one matrices from noisy data, revealing critical thresholds where PCA succeeds or fails in high-dimensional settings.

Contribution

It characterizes the phase transition in eigen-structure for rank-one matrix recovery under Gaussian noise, extending understanding of PCA limitations in high-dimensional regimes.

Findings

Identifies critical noise thresholds for successful recovery

Demonstrates phase transition behavior in eigenvalues

Provides theoretical insights into PCA performance under noise

Abstract

In datasets where the number of parameters is fixed and the number of samples is large, principal component analysis (PCA) is a powerful dimension reduction tool. However, in many contemporary datasets, when the number of parameters is comparable to the sample size, PCA can be misleading. A closely related problem is the following: is it possible to recover a rank-one matrix in the presence of a large amount of noise? In both situations, there is a phase transition in the eigen-structure of the matrix.

Tables

Table 1. Table 1. Relative error averaged over 200 trials. Mean ​ ( E ) Mean 𝐸 \text{Mean}(E) denotes the mean relative error. “Opt” is our optimization procedure and “Eig” is the top eigenvector procedure.

$n$	Opt $\text{Mean}\left(E\right)$	Eig $\text{Mean}\left(E\right)$
500	15.4%	113.5%
1000	14.3%	107.4%
2500	12.5%	123.9%
5000	10.6%	111.8%

Equations

$${\hat{\lambda }}_1\to (1+\sqrt{c}{)}^2\text{mbox}asn\to \infty .$$

$${\hat{\lambda }}_1\to {\lambda }_1(1+\frac{c}{{\lambda }_1-1})\text{mbox}asn\to \infty .$$

$$\underset{n\to \infty }{\lim }|⟨{\mathbf{v}}_1,{\mathbf{x}}_1⟩|=\sqrt{c},$$

$$\underset{n\to \infty }{\lim }|⟨{\mathbf{v}}_1,{\mathbf{x}}_1⟩|=0.$$

$${\hat{\lambda }}_1\to \lambda +\frac{1}{\lambda }$$

$$S_{\mu }(z)={\int }_{\mathbb{R}}\frac{d\mu (t)}{z-t},$$

$${\mu }_{\frac{1}{\sqrt{n}}{M}_n}(x)=\frac{1}{n}\sum _{j=1}^n\delta \left(x-\frac{{\lambda }_j(M_n)}{\sqrt{n}}\right)$$

$$S_n(z)=S_{{\mu }_{\frac{1}{\sqrt{n}}{M}_n}}(z)=\frac{1}{n}Tr\left[{\left(\frac{1}{\sqrt{n}}{M}_n-z{I}_n\right)}^{-1}\right]$$

$$S_{{\mu }_{sc}}(z)={\int }_{\mathbb{R}}\frac{d{\mu }_{sc}(x)}{x-z}=\frac{-z+\sqrt{z^2-4}}{2}.$$

$$\det (z{I}_n-(D+P))=\det (z{I}_n-D)\cdot \det (I_n-(z{I}_n-D{)}^{-1}P).$$

$$Tr(M)=\lambda \sum _{k=1}^n\frac{|u(k){|}^2}{z-{\lambda }_k}.$$

$$\lambda \sum _{k=1}^n\frac{|u(k){|}^2}{z-{\lambda }_k}=1.$$

$${\mu }_n(x)=\sum _{k=1}^n|u(k){|}^2\cdot \delta (x-{\lambda }_k).$$

$$\frac{1}{n}\sum _{k=1}^n\frac{1}{z-{\lambda }_k}=\frac{1}{\lambda }.$$

$$S_{{\mu }_{sc}}(z)=\frac{1}{\lambda }.$$

$$S_{{\mu }_{sc}}^{-1}(\frac{1}{\lambda })=\lambda +\frac{1}{\lambda }.$$

$$0\le \mathbf{x}(i)\le \tau \forall i\in \{1,2,...,N\}$$

$$||\mathbf{x}|{|}_2=1.$$

$${\mathbf{x}}_{\mathbf{k}+1}={\mathbf{x}}_{\mathbf{k}}-\alpha [{\mathbf{x}}_{\mathbf{k}}^{\mathbf{T}}({\mathbf{x}}_{\mathbf{k}}{\mathbf{x}}_{\mathbf{k}}^{\mathbf{T}}-\mathbf{X})+\gamma ({\mathbf{x}}_{\mathbf{k}}^{\mathbf{T}}{\mathbf{x}}_{\mathbf{k}})1^{\mathbf{T}}{]}^{\mathbf{T}}$$

$$\frac{||{\mathbf{x}}_{\mathbf{k}+1}-{\mathbf{x}}_{\mathbf{k}}|{|}_2}{||{\mathbf{x}}_{\mathbf{k}+1}|{|}_2}\le 10^{-5}.$$

$$\text{E}(\mathbf{x})=100\cdot \frac{||{\mathbf{x}}_1-\mathbf{x}|{|}_2}{||{\mathbf{x}}_1|{|}_2}$$

Full text

Keywords: principal component analysis, low rank matrix recovery

1. Introduction

The problem of low-rank matrix recovery has received a lot of attention in the signals processing community over the last 10 years. One reason is the close connection to the matrix completion problem: Given a random subset of the entries of a low-rank matrix, the task is to recover all the entries of this matrix. According to the foundational paper by Candes and Recht [6], this task can be accomplished by solving a convex optimization problem. The matrix completion problem arises in machine learning, image processing, and the Netflix-Prize problem. The practical nature of this problem has motivated many researchers to investigate efficient methods to solve this optimization problem. See, e.g., [5], [11], [12], [13], [14], [10], [8], [18], [7]. This list is by no means exhaustive.

In this paper, we take a different direction. We are interested in the following type of situation: Is it possible to recover a rank-one matrix in the presence of a large amount of noise? A useful data model to keep in mind is the following: $X=\lambda\bf{x}\bf{x}^{T}+G$, where $\bf{x}\in\mathbb{R}^{n}$. The data matrix $X$ represents our observations, and the Gaussian matrix $G$ represents the noise structure. The challenge is to recover the principal vector $\bf{x}$ and the value $\lambda$ from the data matrix $X$. We are especially interested in the asymptotic behaviour of the largest eigenvalue and the leading eigenvector of the data matrix $X$, (as $n\rightarrow\infty$), when the operator norm of $G$ is not negligible compared to the operator norm of $X$. We observe a phase transition phenomenon.

We are also interested in the behaviour of the leading singular vector when the data matrix is a large rectangular matrix, where the number of rows is proportional to the number of columns. Principal component analysis (PCA) is a versatile tool in dimensionality reduction. PCA projects the data onto the principal subspace spanned by the leading eigenvectors of the sample covariance matrix. In theory, these eigenvectors can capture most of the variance in the data. This enables the dimension of the feature space to be reduced, while retaining most of the information. In the contemporary setting, a collection of high-dimensional data can be treated as a low-rank signal with additional noise structure. If the samples of data are organized into a data matrix, then PCA can be used to recover the low-rank signal. It performs well when the number of features, $p$, is small, and the number of samples $n$ is large ([1]). However, in biomedical studies, the number of features $p$ is often comparable to the sample size $n$. In the biomedical setting, the features are measurements on the expression levels of thousands of genes, and $n$ is the thousands of individuals.

1.1. Setting and Motivation

Suppose we have a collection of independently and identically distributed random vectors, $x_{1},x_{2},x_{3},\ldots,x_{n}$ from a $p$-dimensional real Gaussian distribution with mean zero and covariance $\Sigma=\mbox{diag}(\lambda_{1},\lambda_{2},\ldots,\lambda_{M},1,1,\ldots,1)$, where $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{M}>1$. Let $X$ be the $p\times n$ matrix with column vectors $x_{1},x_{2},\ldots,x_{n}\in\mathbb{R}^{p}$. Assume that $0<c<1$ and $\frac{p}{n}=c$. Let $S=\frac{1}{n}XX^{T}$ be the sample covariance matrix. A data scientist wants to know how the largest eigenvalues of the matrix $S$ behave as $n\rightarrow\infty$. Let us consider a specific scenario, when $p=500$ and $n=2000$, is the sample largest eigenvalue $\widehat{\lambda}_{1}$ of the matrix $S$ a good estimator of the true eigenvalue $\lambda_{1}$? That depends on the true value of the largest eigenvalue $\lambda_{1}$.

The following is a simplified version of a theorem of Baik and Silverstein (see [3], [15]):

Theorem 1.1**.**

*Let $\widehat{\lambda}_{1}$ be the largest eigenvalue of $S$.

(1). Suppose $\lambda_{1}\leq 1+\sqrt{c}$ and $\frac{p}{n}\rightarrow c$ as $n\rightarrow\infty$. Then, we have*

[TABLE]

(2). Suppose $\lambda_{1}>1+\sqrt{c}$ and $\frac{p}{n}\rightarrow c$ as $n\rightarrow\infty$. Then, we have

[TABLE]

What these authors observe is that there is a phase transition in the eigen-structure of a matrix when both the rows and columns are large, i.e. when $\frac{p}{n}\rightarrow c,0<c<1$ and $n\rightarrow\infty$.

The phase transition phenomenon can be quite complicated and this has been analyzed in the seminal paper [2]. For other variations on this theme, see, e.g. [16], [9], and [17].

In the data model, $X=\lambda\bf{x_{1}}\bf{x_{1}}^{T}+G$, suppose we have some additional information about the principal vector $\bf{x}$, how can we use that information to recover the vector? We consider the case when each entry of the vector is bounded between [math] and a known constant $\tau$. To be precise, we can take the specific value, $\tau=0.2$. Thus, we know that the vector $\bf{x_{1}}$ lies in a box.

Numerical experiments show that, when $\lambda=4$, the leading eigenvector $\bf{v_{1}}$ of the data matrix $X$ is not a good approximation to the desired vector $\bf{x_{1}}$. In fact, the relative error between $\bf{v_{1}}$ and $\bf{x_{1}}$ often exceeds a hundred percent.

The purpose of this paper is to address both the theoretical and practical aspect of this problem. On the theory side, we observe a phase transition in the largest eigenvalue. Moreover, the result shows that, depending on the true value of $\lambda_{1}$, the leading eigenvector of the data matrix $X$ can be nearly orthogonal to the true vector $\bf{x_{1}}$. This means that some caution is warranted: when $n\rightarrow\infty$, using principal component analysis as an attempt to retrieve the vector $\bf{x_{1}}$ can give a misleading result. In place of a proof to the theorem, we provide a heuristic explanation that gives the main insight to the theorem.

On the practical side, we develop an iterative algorithm to recover the vector $\bf{x_{1}}$ from the matrix $X$. We view this as a box-constrained optimization problem to find a vector ${\bf{x}}$, where the unknown variable satisfies the constraint, $0\leq\|\bf{x}\|_{\infty}\leq\tau$. Compared to the leading eigenvector ${\bf{v_{1}}}$ of the matrix $X$, our algorithm yields a vector that is significantly closer to the desired vector ${\bf{x_{1}}}$.

But first, we need to set some Notations:

$X$ is a symmetric random matrix, $X\in\mathbb{R}^{n\times n}$

${\bf{x}_{1}}$ is a fixed (non-random) vector, $\|{\bf{x_{1}}}\|_{2}=1,$ and ${\bf{x_{1}}}\in\mathbb{R}^{n}$

$G$ is a Gaussian symmetric matrix, $G\in\mathbb{R}^{n\times n}$ and $G=G^{T}$, where

$G(i,j)$ are independent, normally distributed with mean 0 and variance $\frac{1}{n}$ for $i<j$, and $G(i,i)$ is normally distributed with mean 0 and variance $\frac{2}{n}$

The following theorem is the phase transition phenomenon (for symmetric matrices).

Important Note Added: After the initial preparation of an earlier version of this manuscript, we learned that this is a version of a theorem of Florent Benaych-Georges and Raj Rao Nadakuditi, see [4]. We encourage the reader to consult this beautifully written paper. We are grateful to the authors of that paper for bringing this to our attention.

Theorem 1.2**.**

*Let $X=\lambda{\bf{x_{1}}}{\bf{x_{1}}}^{T}+G$, where G is Gaussian symmetric matrix. Pick $\tau=0.2$.

Suppose ${\bf{x_{1}}}$ is a fixed vector of length 1, and $0\leq{\bf{x_{1}}}(j)\leq\tau$, for $1\leq j\leq n$.

Let $\widehat{\lambda}_{1}$ be the largest eigenvalue of the matrix $X$. Let ${\bf{v_{1}}}$ be the leading eigenvector of the matrix $X$, i.e. ${\bf{v_{1}}}$ is the eigenvector that corresponds to $\widehat{\lambda}$.*

Then, if $\lambda>1$, we have

[TABLE]

where $c=1-\frac{1}{\lambda^{2}}$. Otherwise, if $\lambda\leq 1$, we have

[TABLE]

For the largest eigenvalue of the matrix $X$, the following phase transition occurs. If $\lambda\geq 1$, we have

[TABLE]

*as $n\rightarrow\infty$. Otherwise, if $\lambda\leq 1$, we have $\widehat{\lambda}_{1}\rightarrow 2$ as $n\rightarrow\infty$.

Remark 1.3**.**

Interestingly, in their theorem [4], they do not have the hypothesis that the vector ${\bf{x_{1}}}$ satisfied the constraint $0\leq{\bf{x_{1}}}(j)\leq\tau$, for $1\leq j\leq n$, and the conclusion of their theorem remains the same. We include this additional condition in the statement of the theorem, since we explicitly use this condition in our numerical optimization algorithm.

2. Background for Wigner matrices

The symmetric Gaussian random matrix in Theorem 1.2 is an example of a Wigner matrix. We summarize here some background and standard facts regarding the Wigner Semicircular Law for symmetric random matrices. Given any probability measure on the real line, the Stieltjes transform is defined by

[TABLE]

where $z$ is any complex number in the upper half of the complex plane. For any $n\times n$ symmetric matrix $M_{n}$, we can work with the normalized matrix $\frac{1}{\sqrt{n}}M_{n}$ and form its empirical spectral distribution (ESD),

[TABLE]

of $M_{n}$, where $\lambda_{j}(M)$ are the eigenvalues of $M_{n}$. The ESD is a probability measure, also known as the spectral measure for the matrix. For the square matrix $M_{n}$ with spectral measure $\mu(x)=\mu_{M_{n}}(x)$, we can define its corresponding Stieltjes transform. We have the following useful identity,

[TABLE]

where $Tr$ denotes the trace of a matrix, and $I_{n}$ is the $n\times n$ identity matrix. We define the semicircular distribution $\mu_{sc}(x)=\frac{1}{2\pi}\sqrt{4-x^{2}}$. The Wigner semicircular law states that the sequence of ESDs $\mu_{\frac{1}{\sqrt{n}}M_{n}}(x)$ converges almost surely to $\mu_{sc}(x)$. The Stieltjes transform for the spectral measure $\mu_{sc}$ is

[TABLE]

3. Main insight for Theorem 1.2

We now give the heuristic explanation for the quantity $\lambda+\frac{1}{\lambda}$ that appears in equation (1.1) in the phase transition phenomenon of Theorem 1.2.

Note Added: After the initial preparation of an earlier version of this manuscript, we learned that there is a similar discussion in [4], and is accompanied by a rigorous proof. We are grateful to the authors of that paper for bringing this to our attention.

Recall that $X=\lambda{\bf{x_{1}}}{\bf{x_{1}}}^{T}+G$, where G is Gaussian symmetric matrix. Since $G$ is symmetric, we can write $G=U^{T}DU$, where $U$ is an orthogonal matrix and $D=\mbox{diag}(\lambda_{1},\lambda_{2},\ldots,\lambda_{n})$ is a diagonal matrix. Instead of the matrix $X$, we can consider the matrix $UXU^{T}=D+\lambda U{\bf{x_{1}}}{\bf{x_{1}}}^{T}U^{T}$, i.e. a diagonal matrix $D$ plus a rank one positive matrix $P\equiv\lambda U{\bf{x_{1}}}{\bf{x_{1}}}^{T}U^{T}$. The random orthogonal matrix $U$ rotates the fixed vector $\bf{x}$ of length one to a random vector $\bf{u}$ of length one. The intuition is that when $n$ is large, then with high probability, the vector $\bf{u}$ is uniformly distributed on the unit sphere $\{x\in\mathbb{R}^{n}\colon\|x\|_{2}=1\}.$ Hence, each entry $u(k)$ of the unit-length vector $\bf{u}$ is approximately equal to the square root of $1/n$. Fix $z$ and suppose the matrix $(D-zI_{n})$ is invertible. Then, we have the relation,

[TABLE]

Consider the matrix $M\equiv(zI_{n}-D)^{-1}P.$ Then, 1 is an eigenvalue of the matrix $M$ if and only if $z$ is not an eigenvalue of $D$ and $z$ is an eigevalue of $D+P$. Since the matrix $M=(zI_{n}-D)^{-1}\lambda\bf{u}\bf{u}^{T}$ has rank one, so the trace of $M$ is equal to the only nonzero eigenvalue of $M$. On the other hand, we have

[TABLE]

This implies that $z$ is not an eigenvalue of $D$ and $z$ is an eigenvalue of $D+P$ if and only if

[TABLE]

Here, $u(k)$ are the entries of the vector $\bf{u}$. The left hand side of (3.2) is $\lambda S_{\mu_{n}}(z)$, where $\mu_{n}$ represents a weighted spectral measure associated to the diagonal matrix $D$,

[TABLE]

Recall that when $n$ is large, the square of each entry $u(k)$ of the vector $\bf{u}$ is about $1/n$, with high probability. Thus, we replace the previous relation (3.2) with

[TABLE]

But the left hand side of equation (3.3) converges to the Stieltjes transform of the semicircular distribution, so we conclude that

[TABLE]

Inverting the Stieltjes transform, we have $z=S_{\mu_{sc}}^{-1}(\frac{1}{\lambda})$. By direct calculation, we can verify that

[TABLE]

Finally, since $z$ is an eigenvalue of $D+P$, we have shown that $\lambda+\frac{1}{\lambda}$ is indeed an eigenvalue of $D+P$. This completes our heuristic explanation for the appearance of $\lambda+\frac{1}{\lambda}$ in the phase transition of eigenvalues for symmetric matrices.

4. Optimization Algorithm

Our optimization algorithm is based on gradient descent, but includes additional steps to satisfy our problem’s constraints. The input is the symmetric matrix X, which is the observed data with Gaussian noise structure. We use an iterative algorithm to estimate the unknown vector $\bf{x}$, under the constraint (where $\tau$ is an input parameter):

[TABLE]

As the additive Gaussian noise increases the magnitude of the eigenvalues of the observed matrix X, we perform gradient descent on an estimate of the true gradient (i.e. for $\mathbf{x_{1}x_{1}^{T}}$ instead of X) by penalizing the magnitude of the eigenvalues of our recovered matrix: that is, we penalize the trace of $\mathbf{x}\mathbf{x}^{T}$ by adding as a penalty term the $L_{2}$ norm of x. This penalty helps minimize the nuclear norm of $\mathbf{x}\mathbf{x}^{T}$. Our update equation is as follows:

[TABLE]

Above, $\alpha$ is the usual step size, which we experimentally observe to work best when set in the range $[10^{-3},10^{-1}]$, with slightly better recovery results toward the lower end of the range. The second parameter $\gamma$ is the regularization parameter for the $L_{2}$ norm of x, which we set to $10^{-1}$ and do not change. Our gradient descent continues until the following termination condition is met, which is usually satisfied within roughly 50 iterations when $\alpha=0.1$:

[TABLE]

We initialize $\mathbf{x}$ by generating a length $N$ vector of uniform random numbers in $[0,\tau]$ and then dividing it by its $L_{2}$ norm.

This gradient descent procedure, however, does not account for the main constraint in our optimization problem: the box constraint. To satisfy this constraint, after gradient descent reaches the termination condition above, we apply a projection step that mitigates the effect of the additional additive noise in the off-diagonal entries of the observed matrix X. First, we divide $\mathbf{x}$ by its $L_{2}$ norm. Then, we project $\mathbf{x}$ onto the box $[0,\tau]^{n}$ by setting each $\mathbf{x}(i)=\text{min}(\text{max}(\mathbf{x}(i),0),\tau)$. Finally, we again divide $\mathbf{x}$ by its $L_{2}$ norm.

In our experiments, we initialize the true vector $\mathbf{x_{1}}$ by setting a block of 2% of the entries to $1/\sqrt{2(10^{-2})n}$ and dividing it by its $L_{2}$ norm. We keep $\alpha=10^{-1}$. The observed data X is the rank-one matrix $\lambda\bf{x_{1}x_{1}^{T}}$, plus a Gaussian random matrix $G$, as described in the previous section, with $\lambda=4$. We define the relative error as:

[TABLE]

where x is our recovered vector. Regardless of our selection of $\alpha$ in the range above, the standard deviations of the relative error at each $n$ are consistently below $2\%$ for $n\geq 500$. For sizes of $\mathbf{x_{1}}$ ranging from $500$ to $5000$, we observe that the average relative error for the recovered vector using our optimization is substantially lower in comparison to that of using the leading eigenvector. Average relative error is computed over 200 trials (i.e. draws of $G$ and optimization procedures) per $n$. The results of our experiment are displayed in Table 1.