  Non-negative Matrix Factorization (NMF) is a useful method to extract
features from multivariate data, but an important and sometimes neglected
concern is that NMF can result in non-unique solutions. Often, there exist a
Set of Feasible Solutions (SFS), which makes it more difficult to interpret the
factorization. This problem is especially ignored in cancer genomics, where NMF
is used to infer information about the mutational processes present in the
evolution of cancer. In this paper the extent of non-uniqueness is investigated
for two mutational counts data, and a new sampling algorithm, that can find the
SFS, is introduced. Our sampling algorithm is easy to implement and applies to
an arbitrary rank of NMF. This is in contrast to state of the art, where the
NMF rank must be smaller than or equal to four. For lower ranks we show that
our algorithm performs similarly to the polygon inflation algorithm that is
developed in relations to chemometrics. Furthermore, we show how the size of
the SFS can have a high influence on the appearing variability of a solution.
Our sampling algorithm is implemented in an R package \textbf{SFS}
(\url{https://github.com/ragnhildlaursen/SFS}).
