Rank Selection for Non-negative Matrix Factorization

Consider a non-negative data matrix $X$ with $p$ variables (rows) and $n$ observations (columns). NMF ² attempts to approximate the matrix $X$ by the product $TW$ where $T$ is a $p\times k$ non-negative matrix and $W$ is a $k\times n$ non-negative matrix. When $k\ll np/(n+p)$, the dimension is significantly reduced from the original problem. The matrices $T$ and $W$ are called the feature matrix and the weight matrix respectively. Each column of W represents the corresponding observation from $X$ as a non-negative linear combination of the columns of $T$.

Because of the non-negativity constraints, the features and weights selected by NMF tend to be sparse. Furthermore, the features, or types in the matrix T represent parts of the data X. This is different from other dimension reduction methods like principal component analysis or vector quantization, and makes interpretation of results easier. In many contexts, the non-negative combinations can be viewed as mixtures of the types.

NMF has been widely applied in many areas including image data ², automatic music transcription ³, computational biology ^{4 3}, and microbiome data ⁵. The microbiome refers to the whole habitat, including the microorganisms, their genomes and the environmental conditions.⁶ Microbial communities, which consists of a large number of different microorganisms, have been shown to be associated with human health and global nutrient cycling.⁷⁸⁹¹⁰ The non-negative structure of NMF can be neatly interpreted as expressing the microbial community as a combination of subcommunities, which may represent interacting groups of microbes. Jiang¹¹ has applied NMF to marine microbes data to investigate the biogeography of microbial function by extracting a small number of features associated with both microbial protein profiles and sampling location. Ko¹² used NMF to estimate functional microbial interactions from microbiome data.

The approximation $X\sim TW$ can be fitted to minimize either the squared error loss, or more generally the negative log-likelihood of $X$ under some assumed parametric distribution. There are multiplicative algorithms for fitting NMF both for squared error loss and for log-likelihood for Poisson distributed $X$ ¹³. These methods are implemented by the R package NMF, which we use for fitting NMF in this paper. While these methods converge robustly, the solution to NMF is not always unique ¹⁴, and even if it is unique, the likelihood surface can be multimodal ¹⁵, meaning that the optimization can converge to a local optimum. This creates difficulties that we need to overcome in our method.

The only tuning parameter in the NMF algorithm is the rank of the feature matrix, $k$, which is the number of features. The interpretation of $k$ is the number of underlying sub-structures or subcommunities extracted from the data. The selection of $k$ not only affects the performance of NMF but also relates to the interpretation of NMF results. Too small $k$ may lead to key features missing and too large $k$ value may have overfitting issues. It is therefore necessary to select an appropriate $k$ value for the data. The purpose of this paper is to develop a new method for selecting this rank $k$.

The usual approach to NMF rank selection is to use expert knowledge ¹⁶. However, for many analyses in practice, no-one has the insight needed to select the rank. Even if there are experts with a good intuition of what rank is appropriate, being able to support the choice with statistical evidence is valuable.

When the expert knowledge is unavailable, a popular approach for rank selection is to compute some model measures for a range of rank values and choose the rank value according to the measure’s criteria. Cophenetic correlation coefficient and Dispersion are the most commonly used measurements ¹⁷. The Cophenetic correlation coefficient method measures the reproducibility of the clustering of observations based on the weight matrix of NMF for many runs with different initial values. Because of the multimodality of the likelihood or squared errors, runs of NMF with different initial values can produce different results. The idea is that for the appropriate number of features, there should be a clear optimum, so runs with different initial values should give the same solution. If the number of features is too large, then the optimal likelihood or mean squared error will vary between runs, leading to a smaller Cophenetic correlation. The method selects the rank value where the Cophenetic correlation has a steep drop off ⁴. The Dispersion method uses a similar procedure to the Cophenetic coefficient method but calculates Dispersion coefficients instead ¹⁸. These two methods both measure the robustness of the clustering analysis results, instead of the goodness-of-fit of the NMF results to the data, thus could lead to suboptimal solutions ¹⁹. Other ideas include comparing residuals, sparseness and Description Length ²⁰. These methods are all ad hoc in that they select the rank according to the plots based on intuition instead of a solid statistical inference results. The results will be inaccurate when a clear drop off point is not available.

Another popular approach is to use cross-validation. For example, using 2-fold cross-validation to split the data randomly into two halves and apply NMF to the two parts separately. The rank minimizing the average reconstruction error between the two parts is selected for NMF ²¹. Lin¹ introduces a rank selection method based on imputation. They randomly delete $30\%$ of the data entry from the original data each time and treat them as missing values in the NMF procedure. The rank is chosen by minimizing the mean square errors of the reconstruction of the missing entries. Their NMF algorithm to handle missing values is available in the R package NNLM ¹. Compared with the Cophenetic correlation coefficient method, Dispersion method and 2-fold cross validation method, the NNLM imputed mean squared error method estimates the rank more accurately when applied to Normal synthetic data with different properties ¹⁷.

There are also some other rank selection methods such as the Bayesian method ²² and singular value decomposition. The Bayesian method requires an assumption that the rank should be small ²², which is not necessarily true for real data. Singular value decomposition chooses the rank where the singular values become small ¹⁶. The challenge of the method arises when there isn’t a clear drop to zero singular value among all ranks. Moreover, the rank of singular values is not closely related to the rank of non-negative matrix factorization.

In this paper, we develop a goodness-of-fit test based on a deconvolved bootstrap distribution and use the test to choose an appropriate rank for NMF. The application of deconvolution to bootstaps with optimisation error appears to be novel. We then compare our rank selection method with bootstrap distribution without deconvolution and NNLM on a range of simulated data and real data. We find that our method produces accurate and stable estimation of ranks for NMF.

Suppose $L(k)$ is the Likelihood of the factorization using rank $k$ NMF and $L(k+1)$ is the Likelihood of the rank $k+1$ NMF results. We construct a likelihood ratio test. Under null hypothesis $H_{0}$, the test statistic is given by:

Here $\sup$ denotes the supremum and $\ln$ represents natural logarithm.

According to standard statistical theory, under certain conditions, the likelihood ratio statistics asymptotically follow a chi-squared distribution. However, because of the non-negative boundary, the conditions for the asymptotic null chi-squared distributions are not met. Moreover, due to the computational difficulty, it is most common that the global maximization of the likelihood is not achieved by applying the NMF algorithm with a single set of initial values. Thus, the results of a given run of the NMF algorithm are subject to computational measurement errors for the likelihood ratio statistic. For the first difficulty, we can use a bootstrap method to compute a correct null distribution. However, the computational measurement error in the likelihood ratio statistic for the bootstrap samples makes the hypothesis testing method lose its power due to the long tail in the null distribution. We develop a method to get a measurement error deconvolved test to improve the power for the likelihood ratio test in NMF rank selection. This work is also an application of our previous work on a general method for measurement error deconvolution density estimation using penalized likelihood ⁵.

For our hypothesis test, we can use a parametric bootstrap to obtain the null distribution of the test statistic. Under the null hypothesis, the data is a realization of a random sample from a distribution $F$ with mean a linear combination of the $k$ features. The $k$ features and their coefficients are estimated by rank $k$ NMF of the original data. In the bootstrap, samples of random datasets are drawn from $\hat{F}$. These samples have the same dimension as the original data. When the data is Poisson distributed, the mean is the only parameter of the distribution. For Normal data, assuming each entry of the data is independently Normally distributed and shares the same variance, both mean and variance need to be estimated from the original data. The variance can be estimated from the residuals of the $k$ rank NMF model. For example, if the rank $k$ NMF gives $X\approx T_{p\times k}W_{k\times n}$. $TW$ is the estimated mean parameter for the Poisson or Normal distributions. The variance of the Normal distribution can be estimated by $\frac{1}{np}\sum_{i,j}(X-TW)_{ij}^{2}$. The negative entries in the bootstrapped sample are replaced by 0 for Normal distributed data. We apply rank $k$ and rank $k+1$ NMF on each bootstrapped sample to get the test statistic $\lambda^{*}$, from which the empirical null distribution of $\lambda$ can be obtained.

Unfortunately, this hypothesis test with bootstrapped null distribution does not perform well at selecting the rank of NMF since the likelihood surface for NMF is often multimodal, and it is common for the algorithm to converge to a suboptimal factorization. This creates noise in the bootstrapped null distribution, thus reducing the power of the test.

To improve the performance, we can rerun the algorithm with multiple starting values on the original data to choose the best log-likelihood. The mean $\hat{F}$ is estimated as the result of NMF that has the best log-likelihood. For Normal data, the variance is the average of variances estimated from each run of NMF to avoid overfitting. We also run NMF with multiple initial values on each bootstrap sample and choose the best log-likelihoods to build the null distribution. The procedure is shown in Algorithm 1.

The simulation results for this method with different numbers of initial values and different true rank values are shown in Section 3.1.

The major problem with the method described above is that to ensure a good bootstrap distribution, we need to fit each bootstrap sample with a large number of different initial values, which is computationally expensive. We therefore develop a new approach to more efficiently compute the bootstrap distribution. The idea is to treat the convergence error as measurement error in the log-likelihood ratio statistic. We can then use a well-developed deconvolution method to estimate the null distribution of the statistic. Deconvolution is a method to estimate a distribution from data with additive measurement error. In our case, if we run only one set of initial values for NMF for each bootstrap sample for each fixed $k$, then we can directly compute the estimated maximum loglikelihood $l(k)=l_{0}(k)+e(k)$, where $l_{0}(k)$ is the globally maximized log-likelihood value (unobserved) and $e(k)$ is the measurement (convergence) error. Thus

where $\lambda_{0}^{*}$ is the true likelihood ratio statistic without convergence error and $e=-2(e(k)-e(k+1))$. The deconvolution method will estimate the distribution of $\lambda_{0}^{*}$ when $\lambda_{0}^{*}$ is not observable but $\lambda_{0}^{*}+e$ is available.

Most deconvolution approaches are based on the Fourier transformation. While its mathematical solution is very neat, estimation of the characteristic function from data is somewhat unstable, and division by possibly small Fourier coefficients for the error data can greatly magnify errors. So in this paper, we use the P-MLE deconvolution method ²³ which is based on maximizing a penalized log-likelihood of the data (here $\lambda^{*}$) with a smoothness penalty. P-MLE also has the advantage over many other competing methods that its implementation can estimate the measurement error distribution non-parametrically from a pure measurement error sample. This is important because the convergence error is unlikely to follow a standard parametric family of distributions.

Figure 1 shows an example of bootstrapped null distribution of likelihood ratio statistic without deconvolution and after deconvolution. The deconvolved distribution of likelihood ratio statistic has a much smaller vairance than the bootstrapped null distribution of likelihood ratio statistic when NMF with single initial value applied to each bootstrap sample as convergence error is removed by deconvolution.

The sampling procedure for bootstrapped datasets is the same as in Section 2.2. However, for each sampled data, rank $k$ NMF and rank $k+1$ NMF are run only once to get $\lambda^{*}_{0}+e$. In order to get a pure error sample, we perform rank $k$ and rank $k+1$ NMF with multiple initial values for an additional bootstrap sample. Since we are trying to maximize the log-likelihood, we will assume that the largest likelihood from multiple runs of different initial values is the true maximum likelihood. If the number of different initial values is large enough, this should be close to be true. The measurement error samples of $e(k)$ and $e(k+1)$ can be obtained from the largest loglikelihood among different initial values minus the loglikelihoods of all NMF runs. The pure error sample for $e$ can be obtained as negative two times the difference between $e(k)$ and $e(k+1)$ for all different pairs of $e(k)$ and $e(k+1)$. Our method is implemented in the R package for Decon-boot-test, which is available on CRAN with the name Deconboottest at https://CRAN.R-project.org/package=DBNMFrank.

A complete description of the procedure is given in Algorithm 2:

We simulate a range of Poisson NMF scenarios with true rank equal to 2, 4, 6, 8, 10 and 30, and a scenario where there is no NMF structure. For true rank 2, 4, 6, 8 and 10, we simulate 100 observations, while for the true rank 30, we simulate 300 observations. For the no NMF structure scenarios, we simulate 131 observations. For rank 2 and rank 4, we simulate five scenarios to assess the sensitivity of the methods when one of the types is close to a non-negative linear combination of the others. For the higher true ranks, we simulate one scenario for each true rank. For ranks 6, 8 and 10, all features used in the simulation are estimated from the healthy individuals in the Qin dataset²⁴, using Poisson NMF with the given rank. One feature for the true rank 2 simulation and three features for the true rank 4 simulation are also estimated from the healthy individuals in the Qin dataset using NMF with ranks 1 and 3 respectively. The Qin data is a human gut metagenomic dataset extracted from 99 healthy people and 25 IBD patients. For rank 30, the features used in the simulation are calculated by applying rank 30 NMF to Person 2’s gut data in the moving picture data ²⁵. The moving picture data is the most detailed investigation of temporal microbiome variation to date. It consists of a long-term sampling series from two human individuals at four body sites: gut, tongue, right and left palm. After removing rows consisting of all 0’s, the total number of different OTUs is 3131 for Person 2’s gut data. Each feature is normalized to sum to one for all scenarios.

For the rank 4 simulations, the fourth feature is chosen in the space spanned by the other three features and a unit vector with all elements the same, such that the perpendicular from the fourth feature to the plane of the other three features passes through the centroid of these three features. We vary the distance of this fourth type from the plane of the other three to assess the sensitivity of the methods in cases where one type is almost a linear combination of the other three.

For rank 2 simulations, the second feature is on the line in the space spanned by the first feature and a unit vector with all elements the same and in the direction from the first feature to the unit vector.

For rank 2 and rank 4, the further the distance is, the more distinct these features are. The distances are set to be 0, 0.0005 0.0008, 0.001, 0.0015 and 0.002. When the distance is 0, the feature matrix degenerates to a lower rank matrix. Thus, the rank 4 feature matrix degenerates to a rank 3 feature matrix and the rank 2 feature matrix becomes rank 1.

For all Poisson NMF data simulations, columns of the weight matrix are generated independently from a uniform distribution from 0 to 1. We adjust each column of the weight matrix to have a fixed sum sampled from the column sums for the healthy group of the Qin data for true NMF rank 2, 4, 6, 8 and 10. The weight matrix is rescaled to have fixed sum sampled from the column sums for Person 2’s gut data of the moving picture data for true NMF rank 30. For microbiome data, these column sums are referred to as sequencing depth. Multiplying the feature matrix and the weight matrix, we get the parameters for the Poisson distribution. Each entry of the data matrix simulated for the Poisson NMF data is generated from a Poisson distribution with the Poisson mean given by the corresponding entry of $TW$. The estimated results of Boot-test, Decon-boot-test and NNLM in 50 replicates of different scenarios are summarized in Table 1, Table 3 and Table 5. To compare the computational costs of the Boot-test and Decon-boot-test, we recorded the average time required to calculate the estimated rank for 50 replicates in each scenario in Table 2, Table 4 and Table 6. Both methods are coded in R (64-bit 4.0.3) and run on 1 node of the Graham Compute cluster of Canada with 1 Intel Xeon at 2.1Ghz with 2G of RAM and 1 core. The NNLM method is not compared here as it uses a different NMF algorithm while estimating the rank. The efficiency of the NMF algorithm is closely related to the computation cost of the rank selection method.

To compare the ability of Decon-boot-test and NNLM to detect Poisson data without NMF structure, we design a no NMF structure Poisson data simulation. The dimension of the no NMF structure Poisson data is the same as Person 1’s gut data from the moving picture dataset, which has 1864 variables and 131 observations. In each replicate, each row of the mean of the no NMF structure Poisson data is generated from a log-multivariate Normal distribution. The mean of the multivariate Normal is a vector of length 1864 whose entries are independently simulated following a Normal distribution with mean 4 and standard deviation 3. For each sample from the multivariate Normal distribution, the mean is rescaled so that the sum of all means is equal to the log of the total sequencing depth from the corresponding sample in Person 1’s gut data. The covariance matrix of the multivariate Normal is generated by 1864 eigenvalues evenly spaced from $3\times 10^{-7}$ to 1 and orthogonal eigenvectors uniformly distributed over the 1863-dimensional sphere. If the Poisson data doesn’t have NMF structure (which can be thought of as having rank equal to sample size), then rank selection methods should select a high rank (which doesn’t satisfy the NMF rank selection rule $k<\frac{n\times p}{n+p}$  ²). The average estimated ranks and standard deviations for the 50 replicates of the two methods are shown in Table 5.

The Boot-test results in these tables indicate that for Poisson NMF data, the accuracy of Boot-test estimates increases with the number of initial values used for each NMF when $d=0.0005$ and $d=0.0008$ and is stable in most other cases. The computation time of Boot-test is positively related to the number of different initial values used for each NMF ($m$). Decon-boot-test’s computation time is between the time spent by $m=3$ and $m=10$ Boot-test. But Decon-boot-test estimates the ranks more accurately than Boot-test in all scenarios except for true rank 4 with $d=0.0005$. In Table 1 and Table 3, the Decon-boot-test selects the true ranks in more replicates than NNLM especially when the distance from one feature to other features is small, for example $d=0.0005$ and $d=0.0008$. In Table 5, Decon-boot-test is comparable with the NNLM method when the NMF rank is low. For the high rank case, summaries of the estimated ranks show that both methods underestimate the rank, but Decon-boot-test gives a fairly close estimate for the rank, while NNLM estimates a very low rank. When there is no NMF structure, Decon-boot-test selects a high rank, while NNLM selects 1 for all replicates, which is somewhat misleading.

We generate Normal NMF data with true rank equal to 2, 3, 4, 6, 8 and 10. For each scenario, we simulate 50 replicates.

Each element in the feature matrix for the Normal NMF data is independently sampled from a gamma distribution with shape parameter 3 and rate parameter 2, then multiplied by a Bernoulli random variable with parameter $0.7$ to control the sparsity of the feature matrix. All replicates in the same scenario share the same feature matrix. The dimension of the features is 2780 initially. After removing rows with all zero elements in the feature matrix, there are 2544 rows remaining in the type 2 feature matrix and 2721 rows in the feature matrices for other ranks.

The 100 columns of the weight matrix $W$ are generated independently from a Dirichlet distribution with parameters all set to be $1.5$ and rescaled by a factor of 10.

Each entry of the data matrix is simulated following a Normal distribution with mean given by the corresponding entry in the matrix product $TW$ and variance 1. Negative entries in the data are replaced by 0 to generate nonnegative data.

Table 7 shows for Normal NMF data, the boot-test’s accuracy increases when more sets of initial values are used for each NMF. Both the Decon-boot-test and NNLM estimate the true ranks for Normal NMF data more accurately than Boot-test. Decon-boot-test and NNLM are both extremely accurate at estimating the rank for all scenarios. Table 8 shows the average time required by Decon-boot-test to calculate the estimated rank for the 50 replicates in each Normal NMF data scenario is usually between the computational expense of $m=3$ and $m=10$ Boot-test.

The simulation results show that Decon-boot-test is better than Boot-test in estimation accuracy for both Poisson data and Normal data. By using a large number of starting points, Boot-test can achieve similar accuracy to Decon-boot-test, but is far more computationally expensive with no benefit in accuracy. When the true rank is small, the performance of Decon-boot-test for Normal data is comparable with NNML, and for Poisson data, its performance is better than NNML, especially when the features are hard to distinguish from each other. When the rank is large, both methods fail to find the true rank but Decon-boot-test’s estimated rank is much closer to the true rank. When the data is not NMF structured, Decon-boot-test selects very large rank, while NNLM selects rank 1 for all replicates. Thus, Decon-boot-test gives a much clearer indication of the lack of NMF structure, which can be used to better interpret the data.

All three authors contributed to the development of methodology, design of the simulations, interpretation of the real data analysis, and writing of the paper. All three authors contributed to the implementation of the method. All authors read and approved the final manuscript.

The second author is supported by NSERC grant RGPIN-2017-05108. The third author is supported by NSERC grant RGPIN/04945-2014.

The authors declare no potential conflict of interests.