# How to choose an optimal number of latent factors in non-negative matrix factorization?

Given a matrix $\mathbf V^{m \times n}$, Non-negative Matrix Factorization (NMF) finds two non-negative matrices $\mathbf W^{m \times k}$ and $\mathbf H^{k \times n}$ (i.e. with all elements $\ge 0$) to represent the decomposed matrix as:

$$\mathbf V \approx \mathbf W\mathbf H,$$

for example by requiring that non-negative $\mathbf W$ and $\mathbf H$ minimize the reconstruction error $$\|\mathbf V-\mathbf W\mathbf H\|^2.$$

Are there common practices to estimate the number $k$ in NMF? How could, for example, cross validation be used for that purpose?

• I don't have any citations (and actually I did a quick search on google scholar and failed to find any), but I believe that cross-validation should be possible. Commented Aug 11, 2014 at 10:53
• Could you tell me more details about how to perform cross validation for NMF? The K values for the Frobenius Norm will always decrease as the number of K increases. Commented Aug 22, 2014 at 14:01
• What are you doing NMF for? Is it to represent $V$ in lower dimension space (unsupervised) or is it to provide recommendations (supervised). How big is your $V$? Do you need to explain a certain percentage of the variance? You can apply CV after you define your objective metric. I would encourage you to think of the application and finding a metric that makes sense. Commented Oct 29, 2014 at 18:06

To choose an optimal number of latent factors in non-negative matrix factorization, use cross-validation.

As you wrote, the aim of NMF is to find low-dimensional $$\mathbf W$$ and $$\mathbf H$$ with all non-negative elements minimizing reconstruction error $$\|\mathbf V-\mathbf W\mathbf H\|^2$$. Imagine that we leave out one element of $$\mathbf V$$, e.g. $$V_{ab}$$, and perform NMF of the resulting matrix with one missing cell. This means finding $$\mathbf W$$ and $$\mathbf H$$ minimizing reconstruction error over all non-missing cells: $$\sum_ {ij\ne ab} (V_{ij}-[\mathbf W\mathbf H]_{ij})^2.$$

Once this is done, we can predict the left out element $$V_{ab}$$ by computing $$[\mathbf W\mathbf H]_{ab}$$ and calculate the prediction error $$e_{ab}=(V_{ab}-[\mathbf W\mathbf H]_{ab})^2.$$ One can repeat this procedure leaving out all elements $$V_{ab}$$ one at a time, and sum up the prediction errors over all $$a$$ and $$b$$. This will result in an overall PRESS value (predicted residual sum of squares) $$E(k)=\sum_{ab}e_{ab}$$ that will depend on $$k$$. Hopefully function $$E(k)$$ will have a minimum that can be used as an 'optimal' $$k$$.

Note that this can be computationally costly, because the NMF has to be repeated for each left out value, and might also be tricky to program (depending on how easy it is to perform NMF with missing values). In PCA one can get around this by leaving out full rows of $$\mathbf V$$ (which accelerates the computations a lot), see my reply in How to perform cross-validation for PCA to determine the number of principal components?, but this is not possible here.

Of course all the usual principles of cross-validation apply here, so one can leave out many cells at a time (instead of only a single one), and/or repeat the procedure for only some random cells instead of looping over all cells. Both approaches can help accelerating the process.

Edit (Mar 2019): See this very nice illustrated write-up by @AlexWilliams: http://alexhwilliams.info/itsneuronalblog/2018/02/26/crossval. Alex uses https://github.com/kimjingu/nonnegfac-python for NMF with missing values.

To my knowledge, there are two good criteria: 1) the cophenetic correlation coefficient and 2) comparing the residual sum of squares against randomized data for a set of ranks (maybe there is a name for that, but I dont remember)

1. Cophenetic correlation coefficient: You repeat NMF several time per rank and you calculate how similar are the results. In other words, how stable are the identified clusters, given that the initial seed is random. Choose the highest K before the cophenetic coefficient drops.

2. RSS against randomized data For any dimensionality reduction approach, there is always a loss of information compared to your original data (estimated by RSS). Now perform NMF for increasing K and calculate RSS with both your original dataset and a randomized dataset. When comparing RSS in function of K, the RSS decreases with increasing K in the original dataset, but this is less the case for the randomized dataset. By comparing both slopes, there should be an K where they cross. In other words, how much information could you afford to lose (=highest K) before being within the noise.

Hope I was clear enough.

Edit: I have found those articles.

1.Jean-P. Brunet, Pablo Tamayo, Todd R. Golub and Jill P. Mesirov. Metagenes and molecular pattern discovery using matrix factorization. In Proceedings of the National Academy of Sciences of the USA, 101(12): 4164-4169, 2004.

2.Attila Frigyesi and Mattias Hoglund. Non-negative matrix factorization for the analysis of complex gene expression data: identification of clinically relevant tumor subtypes. Cancer Informatics, 6: 275-292, 2008.

• It's not clear why the RSS of random data should be lower than the RSS computed with original data when K is small ? For the rest I understand that RSS of random should decrease more slowly than that on the original Data. Commented Dec 22, 2019 at 14:03
• Re #1: Cophenetic correlation. This is ill-premised and highly subjective, perhaps more so than a scree plot. By measuring "robustness" between random restarts, you are simply measuring the number of discoverable local minima. The larger the rank, the more local minima, generally. Only for very clean data is this a reasonably reliable method. Commented Oct 4, 2021 at 18:41
• Re #2: RSS. I tried using RSS and there is never a cross-over in loss of factorizations between randomized and the original inputs. Yes, it works for spectral decompositions (e.g. singular values), but NMF is fundamentally different and the loss of reconstruction for a randomized matrix will always be greater than a non-random factorization. Commented Oct 4, 2021 at 18:44
• Re #1 + #2: Both of the methods you quote are used for cross-validation of PCA, effectively so. However, NMF is fundamentally different. Thus, the answer given by @amoeba is (likely) the only acceptable method here. Commented Oct 4, 2021 at 18:45

In the NMF factorization, the parameter $k$ (noted $r$ in most literature) is the rank of the approximation of $V$ and is chosen such that $k < \text{min}(m, n)$. The choice of the parameter determines the representation of your data $V$ in an over-complete basis composed of the columns of $W$; the $w_i \text{ , } i = 1, 2, \cdots, k$ . The results is that the ranks of matrices $W$ and $H$ have an upper bound of $k$ and the product $WH$ is a low rank approximation of $V$; also $k$ at most. Hence the choice of $k < \text{min}(m, n)$ should constitute a dimensionality reduction where $V$ can be generated/spanned from the aforementioned basis vectors.

Further details can be found in chapter 6 of this book by S. Theodoridis and K. Koutroumbas.

After minimization of your chosen cost function with respect to $W$ and $H$, the optimal choice of $k$, (chosen empirically by working with different feature sub-spaces) should give $V^*$, an approximation of $V$, with features representative of your initial data matrix $V$.

Working with different feature sub-spaces in the sense that, $k$ the number of columns in $W$, is the number of basis vectors in the NMF sub-space. And empirically working with different values of $k$ is tantamount to working with different dimensionality-reduced feature spaces.

• But the question was about how to choose the optimal $k$! Can you provide any insights about that? Commented Aug 11, 2014 at 9:37
• Your explanation of the NMF factorization makes total sense, but the initial question was specifically about the common practices to estimate k. Now you wrote that one can chose k "empirically" (okay) "by working with different feature sub-spaces". I am not sure I understand what "working with different feature sub-spaces" means, could you expand on that? How should one work with them?? What is the recipe to chose k? This is what the question is about (at least as I understood it). Will be happy to revert my downvote! Commented Aug 11, 2014 at 12:48
• I appreciate your edits, and am very sorry for being so dumb. But let's say I have my data, and I [empirically] try various values of $k$ between 1 and 50. How am I supposed to choose the one which worked the best??? This is how I understand the original question, and I cannot find anything in your reply about that. Please let me know if I missed it, or if you think that the original question was different. Commented Aug 11, 2014 at 13:56
• @amoeba That will depend on your application, data, and what you want to accomplish. Is it just the dimensionality reduction, or source separation, etc ? In audio applications for instance, say source separation, the optimal $k$ would be the one that gives you the best quality when listening to the separated audio sources. The motivation for the choice here will of course be different if you were working with images for instance. Commented Aug 11, 2014 at 14:16
• @MarcoFumagalli by just looking at the similarity/difference between $V$ and $WH$ you will choose a high $k$ because the higher the $k$ the lower the error will be (mentioned in the second answer in the paragraph RSS against randomized data). Commented Aug 14, 2019 at 8:56

In the case of Bayesian NMF (and in general Bayesian MF with the only restriction of linearity of the factor in expectation), we have recently obtained valid closed-formulas when the model is correctly specified (our research also propose some fix to some situations of misspecification).

For a Bayesian NMF using Poisson likelihood, and priors any kind with non-zero mean and variance,

\begin{align} \theta_{ik} &\sim F(\mu_{\theta},\sigma_{\theta}^2), \quad \beta_{jk} \sim F(\mu_{\beta},\sigma_{\beta}^2), \nonumber \\ Y_{ij} & \sim \text{Poisson} \left(\sum_{k = 1}^K \theta_{ik}\beta_{jk}\right), \end{align} we obtain the following formula, which can be calculated with empirical estimates of the variance, expected value, and row-wise, and column-wise correlations $$\rho_1$$ and $$\rho_2$$ (one the results of the paper is that the correlation between two entries has only two possible values depending on them being in the same row or column).

\begin{align} K &= \frac{ \tau \mathbf{V}[Y]-\mathbf{E}[Y]}{\rho_1\rho_2} \left( \frac{\mathbf{E}[Y]}{\mathbf{V}[Y]} \right)^2 \end{align}

and $$\rho = 1 − (\tau_1 + \tau_2)$$.

For more generic Bayesian MF models with any observation model $$F_Y$$ such as

\begin{align} \theta_{ik} &\sim F(\mu_{\theta},\sigma_{\theta}^2), \quad \beta_{jk} \sim F(\mu_{\beta},\sigma_{\beta}^2) \nonumber \\ Y_{ij} & \sim F_Y\left(\sum_{k = 1}^K \theta_{ik}\beta_{jk}\right), \text{ with } \mathbf{E}[Y_{ij}]=\sum_{k = 1}^K \theta_{ik}\beta_{jk}, \end{align}

We obtain a similar result, with an extra dependency on the model specific expectation conditional variance $$\mathbf{E}[\mathbf{V}(Y|\theta, \beta)]$$, namely

\begin{align} K &= \frac{ \tau \mathbf{V}[Y]-\mathbf{E}[\mathbf{V}(Y|\theta, \beta)]}{\rho_1\rho_2} \left( \frac{\mathbf{E}[Y]}{\mathbf{V}[Y]} \right)^2 \end{align}.

One caveat is that that formula is not valid for misspecified models, but it still is a starting point, given that we assume the model has some validity before even seeing the data.

More details can be found in the paper:

Silva, E; Kuśmierczyk, T; Hartmann, M; Klami, A. Prior Specification for Bayesian Matrix Factorization via Prior Predictive Matching. JMLR, 2023.