A Scalable Approach to Estimating the Rank of High-Dimensional Data

Let $X_{n\times p}$ be a matrix with mean-centered columns consisting of a low-rank matrix $S$ plus a random matrix $N$. Without loss of generality, we assume $p>n$ for high-dimensional data. Following the construction of Livan et al. (2018), assume matrix $N$ is formed by $n\times p$ samples from $p$ random variables $[x_{1},...,x_{p}]$ drawn from an unknown joint distribution $p(x)$. Let $N_{ij}$ be independent and identically distributed (i.i.d.) entries with $\E(N_{ij})=0$ and $\E({N^{2}}_{ij})=\sigma^{2}$ (Götze et al., 2004). The scaled inner product of $X$ can be written as

We can specify finding the rank of $X$ in terms of eigenvalue decomposition as follows:

where $u_{X_{i}}$ denotes the $i$th column of the eigenvector $U_{X}$ and $\Lambda_{X}$ is a diagonal matrix with diagonal entries $\{\lambda_{X_{1}},...,\lambda_{X_{n}}\}$ as the non-negative eigenvalues sorted in a non-increasing order. When $n\geq p$, the $\Lambda_{X}=diag\{\lambda_{X_{1}},...,\lambda_{X_{p}},0,...,0\}$, where $n-p$ eigenvalues equals zero. Thus, we can rewrite the scaled eigenvalues of matrix $XX^{T}$ in terms of the following approximation:

There are two challenges in the simplification provided in Equation 3. First, adding $N$ to $S$ can result in the underestimation of the elements in $\Lambda_{S}(S)$ using $\Lambda_{X}(S)$. Second, the eigen vectors of $S$ and $N$ are not generally well-aligned ($U_{N}^{T}U_{S}\not\approx I$). As a result, we provide rigorous justification below and in Appendix A. Given $S$ and $N$ are independent, we have $\frac{1}{p}SN^{T}$ and $\frac{1}{p}NS^{T}$ converge to 0 at a rate of $\frac{\Tr(S)}{p}$ when $p\to\infty$ (see Appendix A). $\Lambda_{S}(S)$ is a rank $K$ diagonal matrix with non-zero entries corresponding to the rank of the signal space. When the signal’s magnitude is large compared to the noise, the corresponding vectors in $U_{X}$ are aligned with $U_{S}$ in the first $K$ dimensions. $U_{X}$ can therefore be regarded as as a perturbation on $SS^{T}$, and $U_{X}^{T}U_{S}$ will be approximately identical in the $K$ leading eigenvectors. The remaining $(n-K)$ dimensions have eigenvalues close to zero. At the same time, $W_{p}=\frac{1}{p}NN^{T}$ is a random matrix converging to $I$ (see Appendix B) and $U_{X}^{T}N^{T}NU_{X}$ can be seen as a rotation on the noise matrix, implying that $U_{X}^{T}U_{N}\approx I$.

When $n\geq p$, we calculate the scaled inner product of $W_{n}=\frac{1}{n}X^{T}X$. Similarly, the scaled eigenvalues of $X^{T}X$ can be rewritten as

The top $p$ eigenvalues of $X^{T}X$ and $XX^{T}$ are the same up to a constant $n$ or $p$ and the rest $n-p$ eigenvalues are equal to zero for $n\geq p$. We only focus on the consistent eigenvalues up to the minimum of $n$ or $p$. The eigenvalues in the null space are not considered in this paper. Similarly, we can demonstrate that $\frac{1}{n}S^{T}N=\frac{1}{n}N^{T}S\approx 0$ and that $\frac{1}{n}N^{T}N$ converges to $I$.

The empirical spectral distribution provides an estimate of the eigenvalues $\Lambda_{N}(N)$ for the $W_{p}$ matrix in high- and low-dimensional cases. We focus on the high-dimensional case where $p>n$ in the in-text derivation and provide a distributional upper bound on the maximum eigenvalue of $\Lambda_{N}(N)$ to determine the magnitude of the signal that can be distinguished from noise. The derivation of the signal-to-noise ratio (SNR) bound for low-dimensional cases can be derived similarly.

The proof of the SNR bound can be summarized as follows: (i) provide the distribution of the diagonal elements of the $W_{p}$ matrix. (ii) Provide the distribution of the off-diagonal elements of the $W_{p}$ matrix and show that the off-diagonal elements converge to zero. (iii) Set the bound for the eigenvalues of the $W_{p}$ matrix by the Gershgorin circle theorem (Gerschgorin, 1931). Intuitively, the eigenvalue of the $j$th column of the $W_{p}$ square matrix lies within a circle of the $j$th diagonal elements at a radius of the sum of the absolute values of off-diagonal elements of the $j$th column (Gerschgorin, 1931). (iv) Further set the bound for the maximum eigenvalue of the $W_{p}$ matrix by the order statistics (Fréchet, 1927).

Assuming that the entries of the noise matrix $N$ are i.i.d. random variables drawn from an unknown joint distribution (without normality assumption) and $p$ is sufficiently large, we first note that by the central limit theorem, the diagonal and off-diagonal entries of $W_{p}$ matrix converge in distribution to normal distributions respectively, as $n$ and $m$ goes to infinity, where $m=n(n-1)/2$ is the number of elements in the upper triangular matrix.

Consider the special case where $N$ follows a standard normal distribution, the variance of the main diagonal elements is $\frac{\gamma^{4}-\sigma^{4}}{p}=\frac{3\sigma^{4}-\sigma^{4}}{p}=\frac{2}{p}$, or similarly by the main diagonal of variance for $W_{p}$ as $\frac{1}{p}(n+1)=\frac{2}{p}$. Thus, the distribution of the main diagonal elements is $\mathcal{N}(1,\frac{2}{p})$, the justification for which is rigorously provided in Appendices B and C. Similarly, the distribution of the off-diagonal elements is $\mathcal{N}(0,\frac{1}{p})$, which is rigorously justified in Appendix D. As proved in Appendix B, the first moment of the $W_{p}$ matrix for the normalized noise matrix is the identity matrix, with the off-diagonal entries converging to zero at a rate of $\frac{1}{\sqrt{p}}$ when $p$ is sufficiently large.

Theorems II.1 and II.3 provide the distribution of the diagonal and off-diagonal entries of the $W_{p}$ matrix. We use the results of order statistics to estimate an upper bound of the maximum diagonal element and the maximum radius of the Gershgorin circle of $W_{p}$ matrix. A similar upper bound can be estimated for the diagonal and the maximum radius of the Gershgorin circle for $W_{n}$ matrix when $n\geq p$.

Equation 3 shows that the rank of the signal can be estimated by first estimating $\Lambda_{N}(N)$ and then finding the zero crossing of the difference $\Lambda_{X}(X)-\Lambda_{N}(N)$ corresponding to where any signal in the signal space cannot be distinguished from noise (see Figure 1).

As a more concrete example, consider the case of plotting scaled eigenvalues of $X$ matrix with 100 samples, 500 features, and a true rank equals to 10 as well as random samples from the MP distribution in Figure 1. The red line plots the scaled eigenvalues of the matrix while the blue line plots the scaled eigenvalues of a similarly distributed noise-only matrix. The goal is to extract all the signal components up to $K$ whose eigenvalues are significantly larger than the expected eigenvalues from the MP distribution. The dimension of the detectable signal is the first eigenvalue index where the corresponding eigenvalue is close to the largest eigenvalue sampled from the MP distribution, meaning it cannot be distinguished from noise. The algorithm is implemented as follows:

An implementation of the dimension, ladle, and Hoff’s method as an open-source package for the R Programming environment can be found at https://github.com/WenlanzZ/dimension. We estimate the corrected variance of the MP distribution by the SNR bound to align the tails of the scaled eigenvalues of $X$ and noise. The test for deviation of scaled eigenvalues can be performed in several different ways. In practice, the univariate Bayesian change point (BCP) analysis, proposed by Barry and Hartigan (1993) has been found to be effective and is implemented in the R package bcp (Erdman et al., 2007). Frequently, there are several change points with their posterior probabilities approximating the maximum posterior probability. In order to guarantee the inclusion of all important signals, we propose a two-step procedure to estimate the dimension $K$ in practice. In the first step, we select candidates for $K$ by calculating the BCP on the posterior probability for a second time (double posterior). In the second step, among those candidates with the highest double posterior, we pick $K$ equal to the largest dimension where its posterior probability is larger than $\delta$ times the maximum posterior probability ($\delta\in\mathbb{R}:\delta\geq 0\text{ and }\delta\leq 1$ determine the range of approximation to the maximum posterior probability; by default set $\delta=0.90$).

One limitation observed in simulation studies when implementing the dimension algorithm is related to R package bcp for calculating the posterior probabilities of the change points. The performance of R package bcp is not always stable at the edges of the sequence. Thus, in some extreme cases, a spike occurs at the edge after a long and flat pattern in the sequence of posterior probabilities. To address this problem, we adopted an alarm system that detects flat and spike patterns and trims off sequence before spikes or after a long and flat pattern.

We compared the performance of the proposed rank estimation procedure to the ladle and Hoff’s method. Consider two signal-plus-noise matrices whose true signal rank $K=3$. For a fair comparison, the first matrix setting is the same as in Luo and Li (2016) and the second matrix setting is similar to Hoff (2007):

We simulated 1,000 independent samples from the $X$ matrix and performed rank estimation with different methods. We allowed each method to run for up to 24 hours for 1,000 simulations. The percentage of correct estimation and elapsed time in seconds are reported in Table 1. For example, in the first case of the $X_{1}$ matrix (sigma separated by commas) with signal rank equal to three and 100 samples and 10 features, dimension had 94.8% of correct estimation with a mean absolute error of 0.056 and an average elapsed time of 0.040 seconds; ladle had 93% of correct estimation with a mean absolute error of 0.396 and an average elapsed time of 0.034 seconds; Hoff’s method had 78.8% of correct estimation with a mean absolute error of 0.367 and an average elapsed time of 0.538 seconds. In the sixth case of the $X_{2}$ matrix (sigma as a single integer) with signal rank equal to three and 10 samples and 100 features, dimension had 28.3% of correct estimation with a mean absolute error of 0.980 and an average elapsed time of 0.041 seconds; ladle had 0% of correct estimation with a mean absolute error of 5.930 and an average elapsed time of 0.044 seconds; Hoff’s method had 0% of correct estimation with a mean absolute error of 7 and an average elapsed time of 1.094 seconds. In Table 1, when we fixed $p$ and let $n$ grow for the $n>p$ setting, the dimension method achieved the highest accuracy and performed stably for both the $X_{1}$ and $X_{2}$ matrices, while the ladle method seemed to work only with its own matrix setting. When we fixed $n$ and let $p$ grow for the $p>n$ setting, dimension exceled in both accuracy and computational efficiency. In terms of computational efficiency, the dimension method is more scalable to large data sets compared to the other methods. In extreme cases where noise overwhelms signal and the signal space is too small to be detectable, dimension underestimated the true rank by a mean absolute error of around 6. When signal is moderate in an extremely high-dimensional setting, dimension achieved an accuracy of 91.1% with a mean absolute error of 0.207 and an average elapsed time of 3.767 seconds while the other two methods are computationally inapplicable to providing an estimation.