Mathematics
Computer Science-2021.gz:143230
Computer Science
8070609ad5df72505de73bb1ef7f3dab95392199
A Scalable Approach to Estimating the Rank of High-Dimensional Data

A key challenge to performing effective analyses of high-dimensional data is finding a signal-rich, low-dimensional representation. For linear subspaces, this is generally performed by decomposing a design matrix (via eigenvalue or singular value decomposition) into orthogonal components, and then retaining those components with sufficient variations. This is equivalent to estimating the rank of the matrix and deciding which components to retain is generally carried out using heuristic or ad-hoc approaches such as plotting the decreasing sequence of the eigenvalues and looking for the"elbow"in the plot. While these approaches have been shown to be effective, a poorly calibrated or misjudged elbow location can result in an overabundance of noise or an under-abundance of signal in the low-dimensional representation, making subsequent modeling difficult. In this article, we propose a latent-space-construction procedure to estimate the rank of the detectable signal space of a matrix by retaining components whose variations are significantly greater than random matrices, of which eigenvalues follow a universal March\u{e}nko-Pastur (MP) distribution.


I. INTRODUCTION AND BACKGROUND
High-dimensional data analysis, in which the number of variables or features in a data set is larger than the number of samples, is common in signal processing, computer vision, and genomic studies (White and Smyth, 2005;Yang, 1996). For these types of data sets, it is generally assumed that there is a low-dimensional and information-rich subspace within the data set (Hoff, 2007). Subsequent methods prescribe a variety of dimension-reduction techniques to find a linear vector subspace that provides a low-dimensional representation of high-dimensional data and retaining those components of the data carrying predictive (or otherwise significant) information. A standard construction for this class of problems considers decomposing an n×p matrix X into a rank K (K p) signal matrix S plus a noise matrix N . The determination of K is paramount for dimension-reduction problems to avoid over-or under-fitting subsequent models (Carreira-Perpinán, 1997). However, due to the challenge of dimension estimation and relatively few effective approaches proposed, many models take rank K as an analyst-specified parameter based on ill-defined or ad-hoc methods.
For example, a standard method of estimating rank K is to retain components that explain 70% of total variance or to plot the decreasing sequence of the matrix's eigenvalues (i.e., scree plot) and having an individual pick the point at which the rate of decrease itself drops significantly (i.e., "elbow") (Yong et al., 2013). The eigenvalue index at which the elbow occurs is denoted as K. An alternative approach (jackstraw) implements a statistical test to identify features with a statistically significant association with any linear combination of principal components (PCs) of interest based on a constructed null distribution by randomly permuting the data set (Chung and Storey, 2015). Although jackstraw is widely applied in dimension estimation, it requires the number of PCs that capture "systematic variation from latent variables" as an input for the algorithm, and the test results from jackstraw only provide a subset of suggested cutoffs from input PCs instead of a clear dimension-estimation cutoff (Chung and Storey, 2015). Furthermore, jackstraw is computationally intensive and time-consuming when applied to large data sets.
Several attempts have been made previously to solve the rank estimation problem. Hoff (2007) proposed using the Bayesian rank-estimation procedure based on the signal-plus-Gaussian-noise model, providing prior and posterior distribution of rank K and posterior estimation of right-and left-singular vectors of matrix X. Posterior inferences can be drawn from expectation E X [p(K|X)] and rank can be determined by the posterior mode of p(K|X). The major weakness that prohibits application of Hoff's method to large data sets (e.g., 1,000 nodes) is the computational complexity of the markov chain monte carlo (MCMC) algorithm (Hoff, 2007). In addition, in the markov chain, Hoff (2007) sampled eigenvalues and the corresponding eigenvectors from the full condition distribution instead of the marginal distribution. As a result, the probabilities of sampling nonzero eigenvalues do not mix well across different ranks of the mean matrix, as eigenvalues and the corresponding eigenvectors are dependent (Hoff, 2007). Recently, Luo and Li (2016) introduced a ladle estimator that combined the decreasing sequence of eigenvalues with the increasing pattern of variability of eigenvector directions. They observed that bootstrap variability of the eigenvectors decrease as the eigenvalues get further apart and increase as the eigenvalues get closer. Luo and Li (2016) proposed the ladle estimator based on this observation. They demonstrated the consistency of the ladle estimator and its superior performance in solving various dimension-reduction problems in a low-dimensional setting. While Hoff's method and ladle have been shown to be effective in several settings, if the focus is on large or highdimensional data, then these two methods may not be the best choice for rank estimation.
We propose a procedure to estimate the rank of the detectable signal space of a highdimensional matrix using the signal-plus-noise model under the assumption that the eigenvalues of the noise matrix follow the universal behavior described in Aparicio et al. (2018).
Under this construction, the expected spectrum of the random matrix N follows a universal MP distribution, which can be compared to the eigenvalues of X (Marčenko and Pastur, 1967). In this context, components carrying detectable signal can be characterized as those whose eigenvalues are significantly larger than expected under the MP distribution. This construction provides a method for estimating the dimension of a data set and extracting those components of the data that may not be attributable to randomness. The major contribution of this chapter is to present this novel and practical method, which results in improved accuracy, robustness, and computational efficiency, especially for high-dimensional data.
The paper is organized as follows: section II describes the framework setup and provides a distributional signal-to-ratio bound. Section III illustrates a procedure of estimating detectable signal dimension and demonstrates its performance in relation to other standard approaches, including the ladle and Hoff's estimators. Section IV presents an application of dimension-estimation method in scRNA-seq data analysis. Section V concludes with a discussion.

A. Signal-Plus-Noise Model
Let X n×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 × 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 ) = σ 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 ith column of the eigenvector U X and Λ X is a diagonal matrix with diagonal entries {λ X 1 , ..., λ Xn } as the non-negative eigenvalues sorted in a non-increasing order. When n ≥ p, the Λ X = diag{λ X 1 , ..., λ Xp , 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 Λ S (S) using Λ X (S). Second, the eigen vectors of S and N are not generally well-aligned (U T N U S ≈ I). As a result, we provide rigorous justification below and in Appendix A. Given S and N are independent, we have 1 p SN T and 1 p N S T converge to 0 at a rate of Tr(S) p when p → ∞ (see Appendix A).
Λ 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 T X 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 = 1 p N N T is a random matrix converging to I (see Appendix B) and U T X N T N U X can be seen as a rotation on the noise matrix, implying that U T X U N ≈ I.
When n ≥ p, we calculate the scaled inner product of W n = 1 n X T X. Similarly, the scaled eigenvalues of X T X can be rewritten as n − p eigenvalues are equal to zero for n ≥ 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 1 n S T N = 1 n N T S ≈ 0 and that 1 n N T N converges to I.

B. Empirical Spectral Distribution
By providing an estimate of Λ N (N ), we estimate the amount of detectable signal that can be distinguished from noise in the region where Λ X (X)−Λ N (N ) ≥ 0. We denote W = N N T and its n non-negative eigenvalues as [λ 1 , λ 2 , ..., λ n ]. The joint density of eigenvalues for the Wishart matrices is as follows: where β = 1 for real entries and σ = 1 2 for standard normal distribution (Götze et al., 2004;Livan et al., 2018). When n ≥ p, we define it as Anti-Wishart ensembleW , where we have n − p eigenvalues equal zero. The joint density function of the Anti-Wishart ensemble is similar to Equation 5, with part of the matrix elements being random and the remaining elements non-random. The determination of these non-random elements are related to the first p rows ofW (Livan et al., 2018). The normalization factor can be calculated by Selberg's integral as where γ = n[α + β 2 (n − 1) + 1] and α = β 2 (p − n + 1) − 1 (Dumitriu and Edelman, 2005;Forrester and Warnaar, 2008;Kumar, 2019;Livan et al., 2018). The empirical spectral distribution of W p is defined as As proven in Götze et al. (2004); Livan et al. (2018); Marčenko and Pastur (1967), the expected spectral distribution EF n (y) and F n (y) converge to the MP distribution in probability with density form: where rectangularity ratio lim n,p→∞ (n/p) → c ∈ (0, ∞). Note that ratio of n/p approximates c and is not fixed. Denote δ(y) as the Dirac delta function and et al., 2004;Livan et al., 2018;Vivo, 2008). When p > n, we have the rectangularity ratio c < 1 and the second term in the density form equals zero according to the index function.
When n ≥ p, we have the rectangularity ratio c ≥ 1 and the second term of the density function does not equal zero only when the eigenvalues equal zero. Note that eigenvalues in the null space are not discussed, as we assume that the signal space is smaller than the minimum of n and p. As a result, the Dirac delta function equals zero for both high-and low-dimensional settings and share the same density form, which is as The empirical spectral distribution provides an estimate of the eigenvalues Λ 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 Λ 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 jth column of the W p square matrix lies within a circle of the jth diagonal elements at a radius of the sum of the absolute values of off-diagonal elements of the jth 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.
Theorem 1. Let N ∈ R n×p with i.i.d. elements having mean 0, variance σ 2 , and the fourth moment γ 4 < ∞. The diagonal entries of the W p matrix ξ j (j = 1, ..., n) converge in distribution to N (σ 2 , γ 4 −σ 4 p ) when p is sufficiently large to have the central limit effect.
Proof. see Appendix C.
Theorem 2. Let N ∈ R n×p with i.i.d. elements having mean 0, variance σ 2 , and the fourth moment γ 4 < ∞. The off-diagonal entries of the W p matrix ε ik (i, k = 1, ..., n and i < k) converge in distribution to N (0, γ 4 +σ 4 4p ) when p is sufficiently large to have the central limit effect.
Proof. see Appendix D.
Consider the special case where N follows a standard normal distribution, the variance of the main diagonal elements is γ 4 −σ 4 p = 3σ 4 −σ 4 p = 2 p , or similarly by the main diagonal of variance for W p as 1 p (n + 1) = 2 p . Thus, the distribution of the main diagonal elements is N (1, 2 p ), the justification for which is rigorously provided in Appendices B and C. Similarly, the distribution of the off-diagonal elements is N (0, 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 1 √ p when p is sufficiently large.
Theorems 1 and 2 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 ≥ p.
Theorem 3. Let ξ j be the main diagonal elements of W p , where j = 1, ..., n. The Fisher-Tippett-Gnedenko theorem shows that for the maximum diagonal element, ξ (p) = max(ξ j ). We have We can set the bound of ξ (n) with the inequality e x ≤ e x 2 + x for x ∈ R such that Proof. see Appendix E.
Theorem 4. Let ε ik be the off-diagonal elements of W p , which are normally distributed with mean 0 and variance γ 4 +σ 4 4p , where i, k = 1, ..., n and i < k. Let R k = i =k |ε ik | be the sum of the absolute values of the off-diagonal elements of the kth column of W p matrix with expectation µ GEV = (n−1) 2 (γ 4 +σ 4 ) 2pπ and variance σ 2 . The Fisher-Tippett-Gnedenko theorem shows that for the maximum entry, The eigenvalue of the kth column of the W p matrix can be bound by the Gershgorin circle theorem. With the distributional upper bound of the ξ k and distributional upper bound of R k , we have the upper bound of the largest eigenvalue of W p matrix, which lies within the Gershgorin discs D(ξ (n) , R (n) ) that are centered at ξ (n) with a radius R (n) of rate of convergence n−1 p for p > n.
Proof. see Appendix F.
elements following a standard normal distribution with mean 0 and variance 1. The diagonal entries of the W p matrix ξ j (j = 1, ..., n) converge in distribution to N (1, 2 p ). The Fisher-Tippett-Gnedenko theorem shows that for the maximum diagonal element, ξ (p) = max(ξ j ). We have we can set the bound of ξ (n) by The off-diagonal entries of the W p matrix ε ik (i, k = 1, ..., n and i < k) converge in distribution to N (0, 1 p ). Let R k = i =k |ε ik | be the sum of the absolute values of the off-diagonal elements of the kth column of the W p matrix with expectation µ GEV = 2(n−1) 2 pπ and variance σ 2 GEV = (n−1)(π−2) pπ . The Fisher-Tippett-Gnedenko theorem shows that for the maximum entry, R (n) = max(R k ). We have where b = α 2 −µ GEV σ GEV . We can set the bound of R (n) with the inequality e x ≤ e x 2 + x for x ∈ R such that Thus, we can set the bound of the largest eigenvalue of W p by D(ξ (n) , R (n) ).

III. ESTIMATING THE DETECTABLE SIGNAL DIMENSION IN PRACTICE
A. Algorithm Equation 3 shows that the rank of the signal can be estimated by first estimating Λ N (N ) and then finding the zero crossing of the difference Λ X (X) − Λ 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. Algorithm 1 Algorithm for rank estimation 1: procedure Algorithm(X n×p ) 2: Calculate eigenvalue decomposition for scaled matrix X to dimension n (n n). Denote eigenvalues as Σ x .

3:
Draw random samples from the MP distribution with a corrected variance. Denote eigenvalues as Σ n .

5:
Calculate posterior probability of change points for deviation (∆ = Σ x − Σ n ) with a prior equals to a decreasing linear sequence from 0.9 to 0.

6:
Set K to equal the largest dimension with the maximum posterior probability.

B. Implementation
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 δ times the maximum posterior probability (δ ∈ R : δ ≥ 0 and δ ≤ 1 determine the range of approximation to the maximum posterior probability; by default set δ = 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.

C. Simulation
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 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.  (Adams et al., 2019). We applied the dimension method to a subgroup of control lung sample (Patient ID 001C) with 2,000 highly variable genes and 127 cells. Figure 2 shows that among all the candidates with maximum double posterior probabilities (see Figure 2b), the eighth dimension is the largest with the maximum posterior probability of change point (see Figure 2a).

V. CONCLUSION
In this chapter, we posed the problem of dimension estimation in a signal-plus-noise regime as the eigenvalue index where the difference between a data matrix and its inherent noise, defined by a random matrix, were close to zero. Based on this regime, we proposed a procedure to estimate the dimension of detectable signal, along with a derivation of distributional bounds for the magnitude of detectable signal required to be distinguished from noise. We implemented the proposed procedure and related methods as per the open-source R package at https://github.com/WenlanzZ/dimension. In simulation experiments, we demonstrated that the accuracy of the dimension method and its speed and computationally efficiency advantages over ladle and Hoff's method, both of which are essential for large and high-dimensional data sets. In high-dimensional cases with extremely low signal, the dimension method provided a robust low-rank estimation while the other two methods provided either zero rank or full rank. Our results showed that the dimension method is more accurate in estimating the dimension while remaining computationally scalable compared to current methods in a wide range of scenarios under both high-and low-dimensional settings.
Besides making performance comparisons with the simulated data, we also demonstrated that community detection in the identified signal subspace was beneficial in cell-type identification when analyzing scRNA-seq data sets. In addition to the applications to genomic data sets, the dimension-estimation method has more general applications, such as digital imaging, speech modelling, and neural networks.
Here we provide detailed proofs of the Equation 3. Base on the central limit theorem, we show that 1 p SN T and 1 p N S T ≈ 0. Given that SN T = N S T , we will focus on the proof for with E(N ij ) = 0 and E(N 2 ij ) = σ 2 . Without loss of generality, we normalize N ij so that it has mean zero and variance one.
Thus we have 1 p SN T ≈ 0 when p → ∞ with a convergence rate Tr S p .
A Scalable Approach to Estimating the Rank of High-Dimensional Data APPENDIX B: Without loss of generality, we standardize the columns of random matrix N n×p (p > n) by subtracting off the column-wise means and standardizing such that σ 2 = 1. We can first derive the first moment of W as follows: where the main diagonals for j = 1, 2, ..., p. and the off-diagonal we have for j = 1, 2, ..., p and i = k.
Thus we have the first moment of the W p matrix as We then derive the second moment of the W p matrix as follows: where the expectation of d th i main diagonal elements is Note that E N 4 id i is the second moment of the chi-square distribution and the off-diagonal elements are with expectation where d i = d j .
Thus we have the second moment of W p as E W p 2 = E 1 p 2 N N T N N T = (1 + 1 p (n + 1))I .
Consider the ith row and kth column of N , for 1 ≤ i, k ≤ p and i < k. We can write the off-diagonal of the elements of the W p matrix at the ith row and kth column as The product of two independent random variables can be written as We rewrite N ij N kj as where a = N ij + N kj and b = N ij − N kj . The moments of a and b is derived as follows: Thus, a and b are independent. So that we can calculate the expectation of N ij N kj as and the variance as Var N ij N kj = 1 16 Var a 2 − b 2 = γ 4 + σ 4 4 .
Thus, for the off-diagonal entry ε ik , the expectation is and the variance as When p sufficiently large to have the Lyapunov central limit effect, we have ε ik N 0, γ 4 + σ 4 4p . (D13)
We can bound the the maximum diagonal element ξ (n) with the inequality e x ≤ e x 2 + x for x ∈ R such that e −e −a ≤ e e −2a − e −a .
Thus, we have E{ξ (n) > α 1 } ≤ e e Theorem 2 indicates that the off-diagonal entries ε ik (i, k = 1, ..., n) and i < k are iid. normally distributed random variables with mean 0 and variance γ 4 +σ 4 4p . Given the distribution of the ith row and kth column of the off-diagonal element of the W p matrix.
We can calculate the sum of the absolute values of the off-diagonal elements of the kth column of the W p matrix as R k = i =k |ε ik |. The absolute value of ε ik follows a halfnormal distribution and the off-diagonal elements are independent, thus we can derive the expectation of the sum of the absolute off-diagonal values as follows: and the variance as Var |ε ik | = (n − 1)(π − 2)(γ 4 + σ 4 ) 4pπ .
The Fisher-Tippett-Gnedenko theorem shows that for the maximum entry of R k R (n) = max(R k ) ∼ GEV(µ GEV , σ GEV , 0) , The cumulative distribution function of the Generalized Extreme Value (GEV) distribution is Φ(α 2 ) = E{R (n) > α 2 } = e −e −b , where b = α 2 −µ GEV σ GEV . We can bound the R (n) with the inequality e x ≤ e x 2 + x for x ∈ R such that e −e −b ≤ e e −2b − e −b .