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

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

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.