Rank Suggestion in Non-negative Matrix Factorization: Residual Sensitivity to Initial Conditions (RSIC)
Authors: Marc A. Tunnell
Zachary DeBruine
Erin Carrier
Venue: Accepted by TMLR ()

Abstract: Determining the appropriate rank in Non-negative Matrix Factorization (NMF) is a critical challenge that often requires extensive parameter tuning and domain-specific knowledge. Traditional methods for rank determination focus on identifying a single optimal rank, which may not capture the complex structure inherent in real-world datasets. In this study, we introduce a novel approach called Residual Sensitivity to Intial Conditions (RSIC) that suggests potentially multiple ranks of interest by analyzing the sensitivity of the relative residuals (e.g., relative reconstruction error) to different initializations. By computing the Mean Coordinatewise Interquartile Range (MCI) of the residuals across multiple random initializations, our method identifies regions where the NMF solutions are less sensitive to initial conditions and potentially more meaningful. We evaluate RSIC on a diverse set of datasets, including single-cell gene expression data, image data, and text data, and compare it against current state-of-the-art rank determination methods. Our experiments demonstrate that RSIC effectively identifies relevant ranks consistent with the underlying structure of the data, outperforming traditional methods in scenarios where they are computationally infeasible or less accurate.  This approach provides a more scalable and generalizable solution for rank determination in NMF that does not rely on domain-specific knowledge or assumptions.

Published in Transactions on Machine Learning Research (04/2025)
Rank Suggestion in Non-negative Matrix Factorization: Residual Sensitivity to Initial Conditions (RSIC)
Marc A. Tunnell
hunnellm@@maaigggunn..dd
Grand Valley State University
Zachary J. DeBruine
debruinz@gusu. edu
Grand Valley State Universily
Erin Carrier
carrieer&gusu..cm edu
Grand Valley State University
Reviewed on OpenReview: https: /t opwneevee..........mm net/ per/formm? id 9x
Abstract
Determining the appropriate rank in Non-negative Matrix Factorization (NMF) is criti- cal challenge that often requires extensive parameter tuming and domain-specifie knowledge. Traditional methods for rank determination focus identilying single optimal rank, which may not capture the complex structure inherent in real-world datasets. In this study, we introduce a novel approach called Residual Sensitivity to Initial Conditions (RSIC) that suggests potentially multiple ranks of interest by analyzing the sensitivity of the relative residuals (e.g., relative reconstruction error) different initializations. By computing the Mean Coordinatewise Interquartile Range (MC1) of the residuals across multiple random initializations, method identifies regions where the NMF solutions are less sensitive to initial conditions and polentially more meaningful. We evaluate RSIC on diverse set of datasets, including single-cell gene expression data, image data, and text data, and compare it against current state-of-the-ai rank determination methods. Our experiments demon- strate that RSIC effectively identifies relevant ranks consistent with the underlying structure of the data, outperorming traditional methods scenarios where they are computation- ally infeasible or less accurate. This approach provides more scalable and generalizable solution for rank determination in NMF that does not rely domain-specifie knowledge or assumptions.
1 Introduction
diweeiinnio models of high-dimensional data foundational for exploratory data analyses. Non- negative Matrix Factorization (NMF) has emerged as one such tool for data decomposition and analysis in various domains, including image processing (Guillamet et al. 2002 Lee & Seung 1999 Liu et al. 2012), text
mining (Hassani et al. 2019 Pauca et al. 2004], and bioinformatics (Devarajan 2008 Gaujoux & Scoighe
2010). By decomposing a non-negative matrix into non-negative factors, NMF often able to extract
mcaningful patterns and components from otherwisc hard-to-interpret datascts (Gillis 2020). However, a
critical challenge in applying NMF is determining the appropriate rank of decomposition (Wang & Zhang 2013), which essentially dictates the number of components to extract from the data.
Traditionally, rank determination methods in NMF have largely focused on identifying a single "optimal" rank using heuristic methods or by leveraging additional knowledge of the distribution of the input data. While often useful, these methods have considerable limitations. They often require arbitrary parameter choices, are sensitive to the initialization, or depend on domain-specific knowledge, which may not always be available or easily interpretable. In this study, we introduce novel approach to rank determination that seeks suggest a number of ranks of interest instead of single optimal rank.

---

Published in Transactions on Machine Learning Research (04/2025)
Multiple ranks can be beneficial as, in most cases, there is single one rank that matters (as NMF is fundamentally not classification algorithm). Typically, the data is simply too large to examine in its raw form, and domain experts are looking for more usable lower rank representations from which they can glean insights. For instance, in bioinformatics, single-cell datasets are large, containing gene expression information. Bioinformatics researchers may glean useful insights from any number of different ranks, but cannot manually look through the full data or all ranks.
Our method, Residual Sensitivity to Initial Conditions (RSIC), is based the observation that the re- construction error of NMF is highly sensitive to the imitial conditions of the factorization. This approach is grounded in multi-rcsolution perspective by considering the stability of a rank's residual to its initial conditions. By doing so, we look to open up new avenues for interpreting NMF results, especially on complex datasets where a single rank may not be able to capture all relevant information.
Our method is designed be general and applicable to a wide range of datascts, without requiring domain- specific knowledge or assumptions. Other methods, such as consensus matrix methods (Brunet et al. 2004
Kim & Park 2007), self-comparison methods (Hubert & Arabie
1985
Grossberger et al. 2018 Sotiras
et al. 2017), and rross-validation-bee
Owen & Perry 2009 Kanagal & Sindhwani 2010 Gilad et al. 2020), have been proposed in the literature to determine the rank of NMF The
approaches (Muzzarelli et al. 2019
vast majority of rank selection techniques have a strong preference for lower ranks, which may not always be appropriate for the data at hand. Our method on the other hand, does not have an algorithmic bias for lower ranks and is designed to suggest multiple ranks. Our methodology is applicable across a wide range of domains and does not rely on domain- speciicc parameters or priori assumptions about the distribution of the data.
This paper is organized as follows. Relevant background information on NMF and the methods against which we compare are in SSection 2 A detailed description of our method is given in Seection 3. A high-level description the datasets we compare on is given in [ecctonn 4 Our experimental sctup is given in Section 5) The results are given in /Section 6 Finally, a discussion and conclusion is given in Secction 7
2 Background
2.1
Non-negative Matrix Factorization
NMF is a low-rank matrix decomposition of an m X n non-negative matrix, A, in which non-negativity constraints imposed in computation of the lower rank factor matrices. Given rank k decomposition, the factor matrices W and H are m X k and k x in dimension, respectively. Although variety of methods
for solving for W and H exist, such as hicrarchical alternating least squares (Kimura et al. 2015 Gillis
& Glineur 2011 or gradient descent (Lee & Seung] 2000), we use Sequential Coordinate Descent (SCD)
Franc et al. 2005 Lin 2007 Hsich & Dhillon 2011] and Multiplicative Update (MU) (Lee & Scung 2000
Lin 2007) exclusively in this study. Additionally, although a variety of objective functions exist, to allow for
easier comparison with other work in this field, we minimize the Euclidean distance between the vectorized representation of A and the vectorized reconstruction, formulated as
2 - H p
where ||-| denotes the Frobenius norm of a matrix. Given this minimization problem, which is subject to non-negativity constraints, and letting Bw A Bи AHT Gw be the gram matrix of W, and G be the gram matrix of H", the SCD update rules can be written in vector form as
H - max H. H., + (Bw).: (H" Gw)ss)
(Gw)i.
for all i 21, ,,,,,, and
Woj max (o.W.,, W., + (Bu). (WGH).
(Gu)ás
2

---

Published in Transactions on Machine Learning Research (04/2025)
for all j E (1, [Lin 2007 Hsich & Dhillon 2011 Franc et al. 2005 Note that the max function
is applied element- and it is assumed that W and H are appropriately initialized prior to computation. Similarly, the MU rules can be written as
H HH B 2 (GwH)) W - W o (B @ (WG#)),
where O and @ denote the Hadamard product and division, respectively (Lee & Seung 2000:
Lin 2007
2.2
Limitations of Non-negative Matrix Factorizations
While NMF has been shown to be useful due to its ability to extract meaningful information, it has short- comings. First, the NMF factorization not unique, making the problem of computing the NMF ill-posed (Gillis 2012). There are typically numerous equally good solutions (e.g. factorizations for which the Frobe- nius norms of the residuals are equally small big). Second, the underlying optimization problem is generally non-convex, meaning there is guarantee that the obtained solution is global minimum (e.g., there is no guarantee that the obtained factorization is the best factorization).
Additionally, it is well understood that NMF is highly sensitive to the initial conditions of the problem, and considerable work has gone into getting around this problem, though it is largely domain-dependen, (Rosales et al. 2016 Yang et al., 2021 Devarajan 2008). Furthermore, as discussed in Langville al. (2014). many NMF initializations, such as the random, centroid, SVD-centroid, etc. exist, all with tradeoffs between performance (both convergence and quality of converged to solution) and computational complexity. Some initializations, such as centroid and SVD-centroid are very costly or rely on having other information from other common, yet expensive, factorizations. More recent work (Esposito) 2021) has looked comprehensively at a wider range of NMF initialization schemne,, grouping them into different categories including random- based, deterministic low-rank approaches, clustering baaeed approaches, and evolutionary-basee approaches. Choice of initialization is still an open question and there is no single accepted best practice.
Fundamentally, NMF is a low-rank decomposition, where the factorization can often be meaningfully inter- preted clustering, Under the view of NMF as algorithm which produces a clustering, the matrix W assigns weights to basis vectors in H, indicating how strongly each data point is associated with different clusters. These clusters, which are represented by the basis vectors of H, often group similar data points together in a meaningful way. However, without additional constraints or sparsity enforcing conditions, the NMF is not a hard clustering algorithm (Kim & Park 2008). Whether used for a soft clustering or not, in order for the factorization to be useful, it is crucial to know which rank decompositions are meaningful. Significant work has been done this effect (as discussed in [Section 2.4 with most studies focusing on classification datascts (for which class labels are known). Typically, this is to ensure that there is "truc rank". While this terminology is commonly used, will instead be more specific and indicate "number of underlying classes". While NMF is often used for clustering (and it is not uncommon to evaluate cluster- ing with classification data), fundamentally the decomposition of basis vectors, even at the same rank as the number of classes, may have no correlation with the classes considered as the true classes and may be meaningless.
2.3 Multi-Rank Rank and Hierarchical Approaches to NMF
In many applications, the notion of a single "best" or "optimal" rank is complicated by the multi-resolution structure that is often embedded in real datasets. This is particularly relevant in single-cell transcriptomics, where an abundance of subtle subpopulations yield nested or overlapping signatures of gene coexpression Kiselev et al. 2019 Wu & Wu 2020). A single fixed rank may fail to capture these hicrarchics, such as Tasic et al. (2016), which advocates for hierarchical decompositions instead of a single static factorization in their transcriptomics analysis. In one notable attempt to resolve complex hierarchies of cell state, (Schwartz et al.) (2020) employed divisive clustering to resolve cell type identity at multiple resolutions, showing recovery of rare cell types that usual resolutions would have been unresolvable from their most similar common homologs.
3

---

Published in Transactions on Machine Learning Research (04/2025)
In Brunct et al. (2004), the NMF concept of "metagenes" was introduced to interpret correlated features. While this was often done at a single rank, they recongnized that these "metagenes" can exist at multiple levels of gramularity, ranging from broad cell-type signatures to fine-grained subclusters. Later studies in single-cell analysis highlight the inadequancy of a single rank for understanding cell mixtures (Kiselev et al., 2019 Wu & Wu 2020). For example, Wu & Wu (2020) argue that discrete cell labels only approximate a true hicrarchical topology, so a low-rank factorization can obscure additional relevant substructures that exist the underlying data. Similarly, (Tasic et al. (2016) utilized iterative clustering to reveal finer subpopulations in single-cell data, which parallels multi-rank factorization approach in NMF.
Methodologically,, multi-rank strategies begin by scanning over a range of candidate ranks (e.g., from 2 to some upper limit) and performing repeated factorizations. One then evaluates each solution for consistency, such as by measuring how often a similar factorization appears under different initial conditions or how stable the factors remain when the data are subsampled. Common metrics include silhouette scores, dispersion of the factor loadings, or residual reconstruction error. For example, Brunct et al. (2004) computes consensus matrices for each rank across multiple runs, and a rank is judged superior if its consesnsus matrix exhibits clear block structure, which indicates consistent and well-separated factors.
Because no single criterion universally captures every nuance of a given dataset, many authors combine multiple diagnostics, such as reconstruction error, stability, and domain-specifie interpretability, to arrive at a suitable rank. KKim & Park (2007) show that even small changes a rank can drastically affect the lactorization's 's interpretability, and they advocate for selecting a rank that yields a balance between parsi- mony and capturing sufficient structure of the underlying data. Meanwhile, Gilad et al. (2020) emphasizes that no single rank determination strategy can be trusted across all datasets. Certain rank determination methods are biased towards low-rank representations of data while other methods may reveal rank that is not particularly relevant to domain scientists.
In practice, this means that domain knowledge is essential when determining whether useful insights can be drawn from a given rank. In the case of single-cell transcriptomics, biologically meaningful metrics, such as whether factors align with known pathways or cell types, may override purely statistical considerations. Consequently, multi-rank scanning, combined with robust evaluation criteria and domain-based validation, offers more comprehensive and useful view than forcing the data into a single, potentially under- over- factored representation.
2.4 Previous Work on Rank Determination
A variety of methods to determine the rank have been proposed in the literature. These can largely be broken into three main categories: consensus-matrix methods, self-comparison methods, and cross-vllidation-bass approaches. A brief overview for each of these methods is provided in this section and further details for the implementations we use for all methods we compare against are given in Section 5
2.4.1 Consensus-Mate Methods
Cophenetic correlation and dispersion coefficient both compute a consensus matrix based on the clustering obtained with NMF, then compute their metric based on this matrix (Kim & Park 2007 Brunet et al., 2004). In both cases, the consensus matrix is computed and then averaged over a number of random starts.
After computing the averaged consensus matrix, both methods compute their metric. Cophenetic correlation then chooses a rank based on when the cophenetic correlation first begins to drop (Brunct et al. 2004). However, it is important to note that what constitutes a drop is dependent on how many ranks are plotted as this affects the range of the y-axis and consequently what appears to be drop. Alternatively, the dispersion cocfficient method chooses a rank based on when the dispersion coefficient is maximized (Kim & Park 2007
In the literature, both methods are often run on a relatively small range ranks around the point at which the authors expect the optimal rank to be. With cophenetic correlation, the optimal rank chosen basod on an extremely small drop in the cocfficient over the range. In fact, one run in Brunct al. (2004) was shown to be inconclusive based the lack of drop on the short range of ranks tested. In our testing,
4

---

Published in Transactions on Machine Learning Research (04/2025)
find the behavior of cophenetic correlation to be erratic all of our datasets when testing higher ranks than was tested in the original implementation. Similarly, we find the dispersion coefficient to be generally increasing as a function of rank after the initial drop. In both cases, these metrics leave the user unsure of what to pick unless choosing to run on a relatively small few ranks. It is not always possible to determine an appropriately small range to check, especially when the rank of the underlying data cannot be determined a priori. Furthermore, the number of ranks included when plotting the cophenetic coefficient changes how stretched the plot is left-to-right between consecutive ranks, which drastically affects how steep a drop appears to be.
2.4.2 Self- Celf-Comparison Methods
Self-comparison methods can be subset into two categories, split validation and permutation comparison. Split validation cuts the input matrix randomly in half, reorders the halves by the similarity of the basis vectors, then computes the similarity between the two halves. This similarity metric can take a few forms but has shown success in the past using Adjusted Rand Index (ART) (Hubert & Arabic 1985 Grossberger jet al. 2018) and inner product. (Sotiras al.) 2017). We choose not compare against inner product it had poor performance real data in the testing performed by Muzzarelli et al (011
Permutation compares the slope of the clbow of the reconstruction error of the factorization of A against the slope of the elbow of the reconstruction error of a permuted version of A. Effectively, this compares the ability of NMF to reconstruct the dataset against the ability of NMF to reconstruct a random matrix of exactly the same magnitude as the original matrix. Effectively, when the slope of the reconstruction error of A is equal to that or greater than the slope of the permuted matrix, no extra information able to be extracted from the original dataset compared to random one.
2.4.3 Non- Non-Caaego Methods
The elbow method is popular technique for rank determination used in cluster analysis. This method involves plotting the residual as a function of & and picking the elbow of the curve as the correct number of clusters to use. The elbow in the graph is where the rate of decrease changes, representing the point at which increasing the number of clusters does not significantly improve the fit of the model. Although the elbow method is at least partially subjective, we compare against it for completeness.
Akaike Information Criterion (AIC) was successful in determining rank on time-series data in Cheung et al. (2015) using modified implementation of NMF. We choose not to compare against AIC as it had poor performance in Gilad et al. (2020) and was otherwise criticized by Ito et al. (2016) for requiring assumptions that do not necossarily hold in NMF. In this study, we focus on a general approach that does not make statistical assumptions regarding NMF or the underlying data.
Similarly, we do not compare against Cai et al. (2022), sequential hypothesis testing method, due to their requirement that the underlying data follows certain distributions.
For similar reasons, we do not compare against the category of Bayesian methods due to their requirements
for a priori knowledge of the distribution for prior estimation (Schmidt et al. 2009 Cemgil 2009).
Additionally, we do not compare against methods utilizing minimum description length (MDL) as they assume a statistical model of the NMF (Yamauchi et al. 2012).
Relevance determination is worthy of mention here but is ultimately not relevant to the discussion in this paper. Essentially, these methods identify relevant clusters given larger rank NMF decomposition (Tan & Févotte 2009).
2.4.4 LuaasVVaiiaaiio-Baaae Approaches
A variety of methods for NMF rank determination using cross- validatii (CV) have been proposed in the literature. These include bi-CV by Owen Perry (2009), and imputation-basod CV by Kanagal Sindhwani (2010p. We choose not compare against bi-CV as the results have been shown practice to be unclear unstable by Kanagal & Sindhwani (2010): Gilad et al. 2020
5

---

Published in Transactions on Machine Learning Research (04/2025)
Imputation-basee based CV is performed by optimizing for W and H given an imputed version of A in which percentage of values are denoted as missing (Kanagal & Sindhwani 2010 The use of a Wold holdout pattern has been shown to be performant and is most widely used (Wold 1978 Kanagal & Sindhwani 2010). We denote withheld values as 1 in the binary masking matrix, M, and otherwise. These values are hidden from computation during optimization. In most implementations, the reconstruction error afterward calculated
as
M о (A - W
(2)
|M 2
In the implementation originally proposed by Kanagal & Sindhwani (2010), this is performed multiple times per rank and then averaged; the rank with lowest mean reconstruction error is chosen. For the purpose of comparison in this study, we will call this method KS-CV, based on the last initials of the authors.
A variety of methods have been proposed that directly build on the work of Kanagal & Sindhwani (2010) such as MADImput (Muzzarell et al., 2019), MSEImput (Muzzarelli et al. 2019), and CV2K (Gilad et al. 2020). Each of these three methods optimizes over every rank of interest a number of times, as was performed by Kanagal & Sindhwani (2010). For MADImput, the Median A Absolute Deviation (MAD) of the rection errors is calculated at each rank and the rank with lowest MAD is chosen (Muzzareli et al. 2019).. We choose not to compare against MSEImput because it performed poorly all but simulated data Muuzzrelli et al., 2019). Differing from the other imput atintation hods described, the authors of Gilad el al.
(2020) compute
the error as the I-omoorm of the error over the masked values, computed against normalized version of the initial matrix which allows them to normalize both W and H, as detailed in Algorithm 2 in their paper. The rank with minimum median reconstruction error calculated as stated is chosen, but is adjusted down based on a correction step determined by a. Wilcoxen rank-sum test GGlad al. 2020
2.5 Complexity
The majority of rank dotermination methods compute NMF using a standard optimization algorithm a number of times, then compule their metric. While the cost of computing these metrics not free, it is typically substantially less than the cost of computing the NMF decomposition. On the other hand, imputation-basee based methods must deal with missing values during the computation of the NMF decomposition, fundamentally altering how the NMF decomposition is computed. The complexity of imputation-basee CV methods is often overlooked but is a significant burden in practice.
Before discussing the complexity of computing an NMF decomposition with missing values, we must first discuss the complexity of a standard NMF decomposition. As is common practice, we only count multipli- cation and division, and we assume naive algorithms for common operations such as matrix multiplication. The amount of work performed at a given rank for a single optimization iteration using SCD when no values are missing is described by
2mnk + 2mk2 + 2nk
(3)
This assumes two Gram matrix computations, one computation of Bw and Bh, and one computation of WG and HT Under the assumption that l & min(m, n), the time complexity is (mmnk). The amount of work required for a single iteration of NMF with MU is slightly more, but results in the same overall time complexity.
The computation of CV with missing values is considerably more complex. The implementation provided by Lin & Boutros (2020), and used by Gilad et al. (2020). implements imputation- CV by creating new Gram matrix for each row column affected by missing values. That is to say that each column, High in the update of requires a different Gram matrix, denoted by G6) and created as
G = Adiag(~M,, W) (diag(~M.,)
didgg((MM.M.))
6

---

Published in Transactions on Machine Learning Research (04/2025)
where denotes element wise logical negation. Similarly, each row, Wa, in the update of W requires different Gram matrix, (( which is created as
G Ma.)) (Hdiag(~M..) H..,
Taking into account the need compute all of these Gram matrices, the missing value computation cost per iteration captured by
2mnk + 2mnk + mk + nk²
(4)
This assumnes the computation of m + n Gram matrices, one computation of Bw and Bи, one computation of W...G
G(I for all (1.2...., m] and (H.,))/G. for all E f1, .n....................... and at least one missing value per row and column. This gives time complexity of (mmmk²), but attention should be given to the two equations which model their behavior. The dominating term of Lquation 3 is contained within Equation 4 meaning that nearly all of the work required to compute the original factorization must be performed addition to the additional gram matrix related work. For datasets of relatively small size, this is not particularly burdensome, but entirely unfeasible on larger datasets as discussed later in SSectioo 6.3. It should be noted that the CV2K implementation from Gilad et al. (2020) has lower time complexity but necessitates the in-memory storage of five additional matrices of equal size to A; this is similarly burdensome for matrices of sufficient size.
In practice, as the computed NMF decomposition is dependent on the init ialization essentially all rank detection methods compute the NMF decomposition for many different initializations. Let a be the desired number of initializations for each rank. Let be the number of iterations. Then, the overall standard NMF decomposition cost for given rank k is
al 2mn2 + 2nkk))
(5)
In contrast, the total CV-based NMF decomposition cost with missing values for given rank le is
(tmm2k + 22mk + 2mk2 + 2mk
3 Residual Sensitivity to Initial Conditions (RSIC)
As previously described in Secction 2.2 NMF is highly sensitive to the initial conditions of the W and H matrices. While many initializations exist, as discussod in [Section 2.2 we limit our scope random initializations in this work. We have found that, even amongst factorizations whose residuals have comparable Frobenius norms, the reconstruction error computed any given factorization. This observation highlights the inherent unpredictability in NMF due to random initializations and can be seen clearly in Figure I This figure shows the delta between the smallest and largest relative reconstruction error at each point in the rank-10 reconstruction of the Faces dataset (described in detail later in Section 4.2.3. Note that the relative reconstruction error for a point (i.j) is computed as
any point may be arbitrarily worse better in
Gi.j Wi.hg + ++..
The maximum delta figure is plotted over the models which had within 10% Frobenius norm of the residual of the median residual model. This figure shows that, even within models that have comparable norms, there exists massive deviation in the reconstruction error when considered coordinatewise. Similar behavior was found to be present in all of the datasets that were tested.
Of interest, these observations reveal a pattern within the variability: at certain ranks, the sensitivity to initial conditions as measured by the deltas in reconstruction error diminishes greatly, forming what can be described as "islands of stability". These ranks, where the variance in reconstruction error is minimal despite different initializations, stand out against nearby ranks that exhibit high sensitivity.
7

---

Published in Transactions on Machine Learning Research (04/2025)
04
Figure 1: The delta between the smallest and largest relative reconstruction error each point in the rank- 10 reconstruction of the 400 faces in the Faces dataset.
We hypothesize that these "islands stability" are not merely random occurrences but could indicative of inherent structure or patterns within the data that are particularly well-captured by NMF at these specific ranks. These stable ranks may correspond with factorizations that more effectively distill the essential features of the data and are less influenced by the initial conditions.
To investigate this, compute number of factorizations at each rank, where the randomly generated initial conditions are related across ranks in scheme described next in [Section 3.1 We then attempt to quantify the Residual Sensitivily to Initial Conditions (RSIC) at a given rank using the metric developed in Section 3.2 Plotted as function of the rank, k we then identify "islands of stability" as potential ranks of interest that should be further investigated. Through general testing, we identify 100 the approximate number of random initializations needed for consistent behavior regardless seed, with more providing no meaningful benefit.

---

Published in Transactions on Machine Learning Research (04/2025)
3.1 Progressive Random Initialization
We implement a progressive random initialization scheme, which allows individual random initializations to be related across ranks, smoothing out variations in reconstruction error across ranks of the NMF factoriza- tion for a single progressive random initialization. For given initialization with maximum rank of interest, kmax: let iii and Hini be random matrices of dimension m X kmax and kmax X n, respectively, whose entries are generated from the uniform distribution on the half interval [0, 1.0). Assume is the desired init be the r-th random initialization of Wmi,, and w(n.) init is a
number of initializations per rank. Let w(r)
copy of the left submatrix (wir) init
Then, for a given initialization, r, and letting kmin be the minimum rank of interest, we have the ordered set tuples of initial matrices, H((rAmmttt1))
Likewise for Hinis let ( be copy of the upper submatrix (H
sk;
(Wi(n Init
(wrttkkaax)) Init
.H (H())
miit
Init
init
Assume set, Sint exists for each progressive random initialization we are performing. For simplicity, we index these sets
(r.,k)
W
(r.k) init
H
(r.k) init
init
(++.k-1) contain copy of the respective Given this indexing, this specifically means that the matrices in Siit
matrices in s(+,k) for all & E {kmin + 1, kmin + 2.... kmax]- This is such that w(*,k-1) contains a copy of the left submatrices of all but the last column of w and H init contains copy of the upper submatrices of all but the last row of (t) Finally, let
init
(r.k) umf
NMF A, s(r.k)
be the resulting factorization of A after performing NMF with a given rank and initialization pairing. This factorization transforms the unoptimized matrices, giving their optimized counterparts. This performed over all progressive random initializations and ranks, giving Smmf-
3.2 Mean Coordinatewise Interquartile Range (MCI)
The goal of this metric is to capture the sensitivity the relative reconstruction error at given rank by measuring the average spread of relative reconstruction errors at each point. This computed by taking the interquartile range (IQR) of the relative reconstruction error each point over the number of initializations.
Let R be matrix of dimension x mn, and assume vec (-) flattens x n matrix to dimension vector in any consistent ordering. The values of R(k) are given by
(k) vec (A
(r(kk))) (r.k) nmf nmf
Z (e + vec (A)),
for all r E {1,2, a} with given rank, k, and where e is the vector of ones. The Mean Coordinatewise IQR (MCI) at each rank, k. is computed as
(C)(2)
mn Σ IQ R
where IQR computes the interquartile range of a vector. This performed for all E (kmin: kiin + k.ax)
Computing the MCI over large number of factorizations requires the storage of all of these factorizations. Although this may be significant amount. of space for large matrices, the total amount of work required is

---

Published in Transactions on Machine Learning Research (04/2025)
significantly less than those methods based cross-validation techniques. Additionally, these matrices may be stored on disk and do not need to be kept in memory.
In this study, we choose ranks that have a relatively low RSIC-MCI value compared to their neighbors of higher rank. To do this, we plot RSIC-MCI as a function of rank and look for ranks in pronounced dips or flat-lining after which the metric increases. To tie break along regions where multiple ranks exhibit the same flat behavior, take the last rank prior to an increase in the metric, as that rank is equally as stable, but further decomposes the data. In practice, may be interested in looking at each of these tied ranks.
3.3 Scalability and Work
The scalability of RSIC is a significant advantage over many existing methods. Like sell-comparison methods, the cost of RSIC is dominated by the computational cost of computing the NMF decomposition. Once an NMF decomposition is computed, the cost of computing the MCI metric is dominated by the cost of computing the residual, which requires one matrix multiplication per initialization per rank. Thus, the total cost to compute MCT metric for a given rank & is (ammk)); and the total cost to compute both the NMF decomposition and the metric for a given rank k is
amnk + at (2mnk + 2mk + 2nk²).
Note that when compared to Equation 5 the added cost of computing the MCT is not the dominant cost, and the overall computational complexity for given rank k remains Oatmmkk). Thus, even including the computational cost of computing the MCI metric, it is still substantial improvement over the cost for imputation-based methods, which require Oatmmnk?) work for a given rank k.
Storage Costs The primary burden of the RSIC method the high storage requirements of storing the factorizations. This storage cost is given as
O (k(m + n))
double precision floating point numbers for given rank k. More precisely, the overall storage cost for all ranks and initializations is given by
Amaax
ar am + n)) =8a + +n) (kimax - kimin + 1) (kmax kimin)
8
kmin
bytes, again assuming double precision floating point numbers. For example, on a real-world dataset we later use: suppose a = 100 on a problem with dimensions m = 11967 and n = 28095. Setting hmin = 2 and kmax = 64, the storage cost is 66.63GB. As the size of the datasets increases, the storage cost will increase linearly with any given dimension, number initializations, or maximum rank of interest. Of course, these matrices can be kept on disk and loaded into memory as needed in later steps of the computation.
Parallelizability A major advantage of RSIC is that it is trivially parallelizable Factorizations at different ranks and with different initializations are computed independently. Although there would appear to be a dependence due to the reliance on the progressive random initialization scheme, it is possible to produce these progressive initializations in any order. Thus, for the a initializations and & ranks, a parallel implementation of RSTC may be performed on ak different machines. Further, the computation of the MCT metric is only partially dependent on the computation of initializations and may begin to be computed as soon as the first initialization complete. For simplicity, we will assume that the computation of the MCI started once all initializations have finished computation.
The computation of the MCI for each rank is independent from other ranks, so we focus only on the computation of the MCI for a single rank. For a given rank, the MCI can be computed in a batched manner, where the metric is computed for a batch of b initializations at a time. These batches can be trivially split over whichever dimension of the original matrix most convenient. We note that is possible (and would potentially lead to better parallel performance) to split the computation of the MCI in 3-dimensional manner the rows, columns, and initializations, but do not explore that in this paper.
10

---

Published in Transactions on Machine Learning Research (04/2025)
For the trivial splitting along columns, and for the first batch of columns ranging from 1 to b, the computation of the residual is given by
Ry.1mb vec A.1.b
(
W k)
(
2 (e vee (A.,1a);;
:,1bb)
for all r E (1.......a). Then computing the partial MCT on the first. batch is given by
MCT( partial
Σ TQR R
The second batch would be processed for ( each batch is completely independent of each other. To combine the partial MCI values, we have
r.mb+ 1:2mb and similar for computing the partial MCI, where
MCI( = Σ (CI partial
for all partial MCI for the given rank k. We note that, when splitting by columns, the vectorization of the matrix must be done column-wise and row-wise when splitting by row.
Thus the computation of the MCI is able to be performed in a parallel batch-wise manner. This adds no overhead in terms of work and requires no commmication between threads except to reduce the final MCI value.
4 Datasets
As discussed in Section 2.2, there is no real concept of a true rank and we instead will indicate the number of underlying classes for each dataset for which it known. Additionally, certain datasets do not have a single number of underlying classes as sub-classes may exist which equally make sense to target the analysis of dataset. We discuss any potential sub-classes in our discussion of each dataset. While presented for comparison with the literature, the number of underlying classes does not necessarily correlate directly, if at all, to the optimal rank of an NMF decomposition. We compare method eight datasets from three different disciplines. We break this section into three subsections, one for each type of dataset used. All datasets are formatted such that each row represents a sample and the columns within row are features. This means that the two single cell datasets, introduced next, are transposed from how they are generally presented the literature.
4.1 Single Cell Datasets
ALL-AML was originally described in Golub et al. (1999) and retrieved using the package provided by
Gaujoux Seoighe (2023). This dataset. is 38 x 5000 dimension, consisting of the 5000 most highly varying human genes in the original datasct and taken from 38 bone marrow samples (Golub al., 1999: Gaujoux & Seoighe 2023). , Of these 38 samples, 27 related to acute lymphoblastic leukemia (ALL) and 11 acute myeloid lenkemia (AML) (Golub et al. 1999). Given the existence of two cancer types, the number of underlying classes is technically 2, though most authors find rank 3 NMF decomposition to be more
informative (Brunet et al. 2004).
The PBMC3K dataset was retrieved from [0x Genomics 20166) and has dimension 2700 x 137144 This dataset consists of 2700 cells and 13714 genes and comes pro-filtered for relevance. This dataset has no consensus factual number of underlying classes but has been shown using domain knowledge to contain 9 distinct cell types Du et al. 2020 which will be treated as the number of underlying classes for our purposes.
11

---

Published in Transactions on Machine Learning Research (04/2025)
4.2
Image-Basee Datasets
4.2.1
Swimmer
The Swimmer dataset was first described in Donoho & Stodden (2003) and retrieved from the package released by (Tan & Févotte (2009 Swimmer has dimension 256 1024 and is synthetic dataset consisting of 256 flattened 32 x 32 black and white stick figure images meant to mimic variety of breaststroke positions. There are total of 16 limb positions, giving a total number of underlying classes of 16.
4.2.2 Full Digits & Dig0246
We retrieved the "Optical Recognition of Handwritten Digits" dataset by Alpaydin & Kaynak (1998) from scikit-learn, a scientifie computing package released by Pedregosa et al. (201)) We use this both as the full dataset and additionally use selected subset for experimentation. We call the full dataset Full Digits and the subset Dig0246, containing only the digits {0, 2, 4, 6}. We subset in this manner for the sake of comparison and consistency with the authors of Muzzarelli et al. (2019), who did not run on the full dataset. We should note that the clustering behavior of NMF appears more clearly separate the classes in Dig0246 than with Full Digits. An example of this clustering behavior will be shown later in Section 6.2 The number of underlying classes are 10 and 4 for Full Digits and Dig0246, respectively.
4.2.3 Faces
We retrieved the "AT&T Laboratories Cambridge Faces" dataset by Cambridge (1994) using the package provided by Pedreeoos et al. (2011). This dataset has dimension 400 x 4096, consisting 400 flattened black and white images of size 64 x 64. The dataset consists of 10 people, with 40 images from each. Given the number of subjects, the number of underlying classes 10. The images were originally 92 x 112 dimension but have been modified in the version retrieved from Pedregosa al. (2011).
4.3 Text-Based Datasets
All text datasets were transformed using bag-of-words model, CountVectorizer, provided by Pedregosa ett al. (01II). The maximum document frequency was set to 95%, minimum document frequency to an integer value of 2, and the stop-words set to english as provided by the implementation. The maximum features was set to 4000 and 30000 for NewsGroup4000 and Web of Science, respectively. All other options were left to their defaults.
4.3.1 NewsGroup4000
We obtainod the NewsGroup-4000 dataset, originally found at Lang (1997), using the packagod provided by Pedregosa al. (201I). This dataset is of dimension 11314 4000 and contains 20 topics ranging from baseball to medical. Many of the 20 topics belong to high-level topics including "comp", "rec", "science", "misc", "talk", "alt", and "soc". In some cases there are additional levels between the high-level topics and the low-level topics. Based on the number of topics, we consider number of underlying classes as 20, although due to the different levels of topics, there is no single number of underlying classes.
4.3.2 Web of Science
We obtained the Web of Science dataset, which was originally described in Kowsari et al. (2017) from Kowsari et al. 2019). The data was preprocessed as previously described using CountVectorizer. This dataset is of dimension 11967 x 28095 and contains 11967 documents from 35 categories and parent categories. Fundamentally, there are two different accurate values for the number of underlying classes in this dataset, 35 and 7.
The package may he found htt.t//wwt......att.e.eeee////aa
12

---

Published in Transactions on Machine Learning Research (04/2025)
5 Experimental Setup
We opt to use only publicly available packages for all significant computations. Due to the large variety of datasets and methods we compare, for computational feasibility, we perform 100 random initializations for each method on each dataset. For each dataset, we let Amin = 2 and kmax = min(m, 64) and perform each method on each rank between kimin and kmax inclusively. Since it has been shown that MU converges to a solution more slowly than SCD (Lin & Boutros 2020j, we use SCD as the optimization routine for all methods in which the option is available. We force all methods using SCD to run to 100 optimization iterations per initialization and rank by setting the tolerance to 1e-16. For the remaining methods in which SCD is not an option, MU is selected and forced to run for 500 optimization iterations due to its slower convergence. For consistency with Muzzarelli et al. (2019), all methods optimize the Frobenius norm defined in Equation 1 As mentioned previously, we limit our comparison to subset of metrics and initializations initializations for the reasons described in [Section 2. Before the start of computation for each method and dataset, the random state is set to the seed 123455789. For the computation performed with the packages provided by Gilad et al. (2020) and Lin & Boutros (2020), which offer multi-threading support, computation is performed on a dedicated workstation with a 32 thread Threadripper processor and 64 GB of memory. All other computation is performed sequentially on a consumer-grad, Intel processor with 16 threads and 128 GB of memory.
5.1 Elbow
For the comparison against the clbow method, we determinc the elbow based on the average of the recon- struction error at each rank. The results are provided based on a visual determination of where the bend appears to be. We utilize averaging here to smooth the plotted curve, which makes determining the elbow reasonable task. In our testing, we find that there is little difference between using the mcan or median.
5.2 Cophenetic Correlation & Dispersion Coefficient
For the comparison with cophenetic correlation and dispersion coefficient methods (Brunet et al. 2004 Kim & Park 2007), use the NIMFA package provided by Zitnik & Zupan (2012). This package does not have support for SCD, so use the MU optimizer instead. The initialization type is set to random, the number of initializations to 100, and the number of iterations to 500. The update type is set to euclidean and the objective function to fro. We pass in the range of ranks from kmin to kmax and plot the cophenetic correlation coefficient and dispersion coefficient as determined by the package. We then choose the optimal rank for each based on the critoria set forth by the original authors in Brunet et al. (2004) and Kim & Park (2077) and briefly described in this paper in Section 2.4.1
5.3 Permutation
For the comparison with permutation, we perform the optimization of NMF using the implementation provided by Pedregosa et al.] (201T). We compute the Frobenius norm of both the permuted and non- permuted matrix as a function of rank. Using SciPy's implementation of the Savitzky-Golay filter (Virtanen et al., 2020 we approximate the slope of the residuals as function of rank. We look for the point where the reconstruction error of the unpermuted matrix decreases less sharply than that of the permuted matrix. This corresponds to the point at which the slope of the residuals of the non-permuted matrix is greater than or equal that of the slope of the residuals of the permuted matrix because both are decreasing and the slopes are negative. In terms of the approximated derivative, the rank selected is the point immediately before overtaking the elbow of the permuted matrix, or the point which they are exactly equal. Due to floating point arithmetic, equality is considered a relative tolerance of x 10- computed relative to the larger of two derivative estimates. This is performed 100 times and the median result is computed.
13

---

Published in Transactions on Machine Learning Research (04/2025)
5.4 ARI Method
For the comparison against the ART method, we randomly split. the matrix into two equal parts along the dimension relating to factors. If splitting into even parts is not possible, the larger split is truncated to size. Then NMF is computed on cach split using the implementation provided by Podregoss et al. (2011). The cosine distance between factors is computed using the package provided by Virtanen et al. (2020), giving a distance matrix. This is then passed to the linear_sune assignment function as provided by Virtanen itt al. (2020), which implements the Hungarian algorithm described in Crouse (2016). The result is the ARI, which is then averaged over all initializations a given rank. The rank with highest mean ARI is selected. We note that ran with both the mean and median ARI, which resulted in the same decided rank in all cases.
5.5 Cross Validation Methods
For the comparison against KS-CV and MADImput, use the package provided by Lin & Boutros (0220)
The output of this package is the reconstruction error at each rank as previously defined in Equation 2 Using this output, we are able to compute each of these metrics and choose a rank based on where it is minimized. For the comparison against CV2K, we use the package provided by Gilad et al. 2020 which returns both a similar output and the chosen rank based on their criteria.
We note that we have chosen use a different initialization scheme than was used by Gilad al. 2020 and Muzzarelli et al. (2019). In Muzzarelli et al. (201)), the authors performed 100 runs, initializing 20 times in each run before choosing the best model as defined by reconstruction error. In Gilad et al. (2020), the authors initialize large number of times, optimizing for hundreds of iterations, and chose the best model defined by reconstruction error to continue. In order to perform this, we would need to run 2000 times for the Muzzarolli et al. (2019) method and 10000 times for the Gilad et al. (2020) method In addition to the considerable length of time required, we do not believe it is fair to arbitrarily run one method an order of magnitude more times than another. Instead, we focus on running all methods for the same number of initializations. For methods that rely differences across random initializations, this means that they do not get the added benefit of testing against thousands of initializations.
In order to enforce these requirements, tolerance may be provided to the package by Lin & Boutros (020) For the package by Gilad et al. (2020), minor modifications to the code are required. The code is modified to force running all 100 itcrations and the intt_facter catt_factormmttricee matrices function is modified to return the first generated W and H.
6 Results
In this section, we present results and compare against other methods. In our analysis of the PBMC3K, NewsGroup4000, and Web of Science datasets, significant computational hurdles were encountered. Due system, we are forced to run CV2k with only 16 threads for both the PBMC3k and NewsCroup44000 datasets. This limitation results in an excessively prolonged runtime of over 30 days for PBMC3k and longer for NewsGroup4000. Given the size of the Web of Science dataset, a further reduction in threads would be necessary, leading to an even longer computational burden.
RAM constraints
NNLM-CV not constrained by memory requirements but is unable to complete even a single run over all ranks on the PBMC3k dataset in 24 hours, meaning the full 100 runs would take over 100 days. The larger scales of NewsGroup4000 and Web of Science implies that completion times would be significantly longer. This difficulty, owing to the incredible time complexity of the method as described in Equation 4 lead us to the conclusion that the burden of computation for these methods is infeasible.
Consequently, the result table for these methods on the mentioned datasets are denoted as "N/A". " In other instances, where the output was inconclusive, we have designated the results "undetermined" (Und.). For instance, undetermined is noted when using the elbow method and there is no clear elbow or when using
2 A number of initialization was not described in the CV2K paper but was instead found on the Github linked in the paper.
14

---

Published in Transactions on Machine Learning Research (04/2025)
the permutation method but the slope of the permuted is steeper than the slope of the unpermated data from the beginning. The remaining results are presented either integer values or in increasing order set format, which is applicable only to permutation as well as MCI-RSIC as previously described.
6.1 Single Cell Datasets
The resulting suggested ranks for cach of the single-cell datasots are shown in Table 1 MCI metric as & function of rank on the ALI-AML dataset in Figure 2 The horizontal axis, representing rank, ranges from 2 to 38. Based on our selection criteria described before, we find rank 5 to be the first island stability, Note that ranks through all have the same or similar MCI values, but rank 5 is the first rank before significant increase in MCI.
We show the RSIC-
We find that method and permutation perform poorly the ALI.AML. dataset, whereas cophenetic, dispersion, and CV2k all converge the generally agreed upon rank and ARI and KS-CV converge the number of underlying classes. For PBMC3K, no method returned a solution that contended with the number of underlying classes, although MCT-RSIC and elbow both return nearby ranks.
'Recall this dataset has no factual underlying classes.
0.5
Mean Coordinatewise IQR (Relative)
23
0.4
0.3
19
0.2
30
10
0.1
0.0
2
17
23
28
Rank of Factorization
Figure 2: Mean Coordinatewise IQR (MCI) (y-axis) vs rank (æ-axis) for the ALL-AML dataset for ranks 2 through 38. We identify ranks 5, 10, 19, 23, and 30 as "islands of stability" and thereby potential ranks of interest.
Table 1: Results for all evaluated rank determination methods in comparison to the number of underlying classes the ALL-AMI. and PBMC3K single cell datasets.
Method
ALL-AML
PBMC3k
MCI-RSIC
{5,
10, 19,
23,
30)
{3. 7,
11}
Elbow
Cophenetic
Dispersion
Permutation
41
ARI
KS-CV
N/A
CV2K
N/A
MADImput
N/A
Number of Underlying Classes
g3
15

---

Published in Transactions on Machine Learning Research (04/2025)
6.2 Image-Basees Datasets
We show the RSIC-MCI metric function of rank the Full Digits dataset in Figure This image plots the MCI metric from rank 2 to 64. Like the other methods tested, this method performs equivalently poorly on the Full Digits dataset when compared to the number of underlying classes. There is significant mixing between the ranks, which is indicative poor clustering behavior.
In addition to the RSIC-MCI metric on Full Digits, we show the metric on Dig0246 in Figure 4 This image plots the RSIC-MCI metric from rank 2 to 64. This clearly shows island of stability rank and provides evidence for one at rank 6.
0.35
Mean Coordinatewise IOR (Relative)
0.30
22
0.25
0.20
0.15
0.10
3
0.05
0.00
2
11
20
29
37
46
55
64
Rank of Factorization
Figure 3: Mean Coordinatewise IQR (MCI) (y-axis) vs rank (I-axis) for the Full Digits dataset for ranks 2 through 64. We identify ranks 11, and 22 as "islands of stability" and thereby potential ranks of interest.
0.30
Mean Coordinatewise IQR (Relative)
0.25
0.20
0.15
0.10
0.05
4
0.00
11
29
37
55
Rank of Factorization
Figure 4: Mean Coordinatewise TQR (MCI) (j-axis) rank (z-axis) for the Dig0246 dataset. for ranks 2 through 64. We identify ranks 4 and 6 as "islands of stability" and thereby potential ranks of interest.
16

---

Published in Transactions on Machine Learning Research (04/2025)
0.200
Mean Coordinatewise IQR (Relative)
0.175
0.150
0.125
0.100
16
0.075
0.050
0.025
2
11
20
29
37
46
55
64
Rank of Factorization
Figure 5: Mean Coordinatewise 1QR (MCI) (y-axis) vs rank (æ-axis) for the Swimmer dataset for ranks 2 through 64. We identify rank 16 as the only "island of stability" and thereby potential rank of interest.
(a) Rank 10
(b) Rank 11
Figure 6: Clustering behavior of the Full Digits dataset at ranks 10 and 11.
Finally, we show the RSIC-MCT metric as a function of rank on the Swimmer dataset in Figure 5) This image plots the MCI metric from rank to and clearly shows an island of stability at rank 16, which corresponds to the number of underlying classes in the dataset.
Figure 6 shows poor clustering behavior at the rank of 10 (number of underlying classes) the Full Digits dataset. Additionally, show the clustering behavior at rank 11, which is nearby rank that we suggest. This figure shows that the clustering behavior at both ranks is rather poor, and provides evidence for why all methods perform poorly on this dataset. Indeed, there appears to be no clear separation between the classes the dataset. This poor clustering behavior is present in each tested rank, and is indicative of the poor performance of all methods on this dataset.
For the image based datasets, MCI-RSIC and elbow return the true number of underlying classes on the Swimmer dataset. The permutation method consistently detects rank greater than the number of un- derlying classes, across all image datasets. For the faces dataset, MCI-RSIC identifies ranks fundamentally different than the number of underlying classes and only elbow chooses a rank matching the number of un- derlying classes. On Full Digits, all methods perform poorly though the clustering behavior of this dataset is poor as seen in Figure 6 For Dig0246, MCI-RSIC, elbow, cophenetic, dispersion, and ARI all identify a rank matching the number of underlying classes. We note that MADImput found the underlying number of classes their paper and the discrepancy here is most likely due to their use of 1900 more initializations than allowed in our study (Muzzarelli et al. 2019h. The results shown in Table 2
17

---

Published in Transactions on Machine Learning Research (04/2025)
Table 2: Results for all evaluated rank determination methods in comparison to number of underlying classes on Swimmer, Faces, Full Digits, and Dig0246 image datasets.
Method
Swimmer
Faces
Full Digits
Dig0246
MCI-RSIC
16
(2,
6,
9, 18)
f3, 11, 22)
(4, 6)
Elbow
16
(Und.)
Cophenetic
Und.
2
Dispersion
64+
2
Permutation
18
50
22
17
ART
13
2
5
4
KS-CV
14
51
12
12
CV2K
17
16
16
MADImpul
13
2
Number of Underlying Classes
16
10
10
Table 3: Results for all evaluated rank determination methods in comparison to the number of underlying classes on NewsGroup4400 and Web of Science text datasets.
Method
News Group40000
Web of Science
MCI-RSIC
{4, 12, 24, 463
{3, 7,
11,
20, 39)
Elbow
9
6
Cophenetic
5
3
Dispersion
12
3
Permutation
Und.
64+
ARI
3
KS-CV
N/A
N/A
CV2K
N/A
N/A
MADImput.
N/A
N/A
Number of Underlying Classes
20
[35,7] 7
6.3 Text-Based Datasets
No method selects rank matching the true number of underlying classes for NewsGroup4000, although MCI-RSIC and dispersion both return nearby rank. The others underestimate more substantially (elbow, cophenetic, dispersion, ARI), are infcasible to run (KS-CV, CV2K, MADImput), or failed to identify any rank (permutation).
For Web of Science, MCI-RSIC identifies the number of underlying classes for the category as rank of interest and nearby rank to the number of subcategories. Elbow detects nearby rank for the category. All other methods either substantially undershoot (cophenetic, dipsersion, ARI), are infeasible to run (KS- CV, CV2K, MADImput), or indicated the rank was greater than the tested range (permutation). These results are shown in TTble 3
7 Discussion & Conclusion
In this paper, introduced RSIC, novel method for determining ranks of interest in NMF, Unlike tra- ditional methods which aim to identify single optimal rank-ofter requiring extensive parameter tuning and domain-specific knowledge-our approach identifies multiple, possibly relevant ranks by analyzing the sensitivity of the reconstruction residual to different initial conditions (random initializations in the case of this paper). This allows for more nuanced understanding of the data's underlying structure and provides some flexibility exploratory data analysis.
Our method identifies "islands of stability", which are ranks where the NMF solutions are less sensitive initialization and, therefore, more likely to represent meaningful decompositions of the data. We quantified
18

---

Published in Transactions on Machine Learning Research (04/2025)
this stability using the Mean Coordinatewise Interquartile Range (MCI) of the relative reconstruction error across multiple initializations. By doing so, we highlighted ranks that consistently produce stable and interpretable factors, providing insights that single-rank methods may overlook.
We evaluated RSIC on a diverse set of datasets across various domains, including single-cell gene expression data, image datasets, and text corpora. Our experiments demonstrated that RSIC effectively identifies ranks of interest that are consistent or close to the number of underlying classes in the data. Of note, our method performed well on large-scale datasets where other methods tended to undershoot or where cross- voss-validation- based approaches were infeasible due to their high computational complexity.
Comparative analysis with existing methods, including consensus-matrix methods like cophenetic correlation coelficient and dispersion coelficient, sell-comparison mots like the ARI, and cross-validation approaches, showed that RSIC is competitive and often superior in identifying meaningful ranks.
However, our method is not without limitations. In datasets where the underlying structure less pro- nounced or when the data does not exhibit clear stability islands, RSIC may suggest multiple ranks, re- quiring further analysis to select the most appropriate one. For example, in the ALL-AML dataset, RSIC showed perfect stability for all ranks at or below 5, but we select 5 based on our selection criteria-this is a limitation of our method. Additionally, while our approach reduces the computational burden compared to some methods, it still requires multiple NMF computations across a range of ranks and initializations.
For future work, we plan to refine the RSIC metric to better handle datasets with subtle or hierarchical structures. Incorporating additional criteria, such as sparsity constraints or domain-specifte knowledge, may help to further refine the ranks suggested by RSIC. Additionally, other types of initialization schemes, such those based clustering, dimensionality reduction techniques, or other randomization schemes, to further explore the space of possible initializations. Further, the method could benefit from a window-based smoothing (e.g., Savitzky-Golay of the MCI values to reduce the noise in the output, which could be particularly uscful in datasets with high variability in the reconstruction error. We also aim to explore the theoretical underpinnings of the observed islands of stability to provide deeper insights into why certain ranks yield more stable decompositions.
believe that RSIC could benefit from
RSIC offers robust and significantly more scalable approach for rank suggestion in NMF, taking steps toward bridging the gap between the need for meaningful data decompositions and the practical constraints of computational resources. By providing a selection of relevant ranks and highlighting areas of stabil- ity, our method empowers practitioners to make more informed decision in exploratory data analysis and dimensionality reduction tasks.
Broader Impact Statement
The proposed RSIC method enhances the applicability of NMF in various fields by provide a more reliable and computationally efficient means of determining relevant ranks. This can benefit areas like bioinformatics, image processing, and text mining, where understanding the underlying data structure is crucial. However, with any data analysis tool, care must be taken into considering when drawing interpretations drawn from NMF decompositions are valid and do not inadvertently reinforce biases present in the data or researchers.
Author Contributions
Both Mare Tunnell and Erin Carrier jointly contributed to the conceptualization of the project and different methods used in the rank suggestion method and both were involved in conceplualization of the experiments and analysis of the results. Mare Tunnell implemented the method, implemented the experi- ments, and drafted the manuscripts. Zachary J. DeBruine provided insight into NMF and valddaaaioo and was crucial in funding acquisition. Erin Carrier provided guidance on the methodological approach, oversaw the project, edited the paper for clarity and coherence.
19

---

Published in Transactions on Machine Learning Research (04/2025)
Acknowledgments
This project was in part funded by the Chan Zuckerberg Initiative Data Insights Grant 24F2022-249404. This research made use of computational resources the Distributed Execution Network (DEN) lab at GVSU.
References
10x Genomics. Peripheral blood mononuclear cells (PBMCs) from a healthy donor (same donor as pbmc6k).. July 2016.
URL cf. cf.10xgenomics cs. wttpaiic.ccoceeeceeeeeeeeeeeeeee http://aappassellllp////////aa/ http:p/mmm...ppoo......eeeeeee
filtered _enne __c bc_ sismattrices..
tar.gz gz
E. Alpaydin and C. Kaynak. Optical recognition of handwritten digits, 1998. URL https://dbi.org/// 24432/C50P49
Jean-Philippe Brunet, Pablo Tamayo, Todd R. Golub, and Jill P. Mesirov. Metagenes and molecular pattern discovery using matrix factorization. Proceedings of the National Academy of Sciences, (01(4)):414 4169, March 2004. doi: 10.1073/ 00.1773/paa.0008850 0308533005..1 URL httpp:://ao..orgg////// 10t.://.00///11111/ppppp....m http://ppaaa.00531100000 0308531101
Yun Cai, Hong Gu, and Toby Kenney. Rank selection for non-negative matrix factorization. 2022. URL
https://api. stw.:piiissmmaanncchhchoaaaa www.oooppg............. oorpaartcsccccoooooooooooo/ooooo 253254815
AT&T Laboratories Cambridge. The database of faces, April 1994. URL https://cam-orl.com uk/
facedatabase. html]
Ali Taylan Cemgil. Bayesian inference for nonnegative matrix factorisation models. Computational Intelli- URL
gence and Neuroscience, 2009:1-17, 2009. doi:
2009/785152
Vincent C. K. Cheung, Karthik Devarajan, Giacomo Severini, Andrea Turolla, and Paolo Bonato. De- composing time series data by non-negative matrix factorization algorithm with temporally con- strained coefficients. In 2015 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEFE, August 2015. doi: 10.1109/embc..215.733333 1010999embc 2015. 7319146
URL https:
David Crouse. On implementing 2D rectangular assignment algorithms. IEEE Transactions on Aerospace and Electroníc Systems, 52(())1667- 2016. doi: 10.1109 TTw.TT..................
Karthik Devarajan. Nonnegative matrix factorization: An analytical and interpretive tool in computational biology. PLoS Compulational Biology, 4(7):e1000022) July 2008. ISSN 1553-7358. doi: 10.1371 10.1371/journal /journal. 1000.1000029 URL http::/dxx... doi orgo 1tt.:/00..////ouuuruu http://///oooaa......ooo........ pcbi000009pp 1000029
David Donoho and Victoria Stodden. When does non-negative matrix factorization give a correct decomposi- tion into parts? In S. Thrun, L. Saul, and B. Schölkopf (eds.), Advances in Neural Information Processing
Systems, volume 16. MIT Press, 2003. URL htttps:/prooceeing..et www.ewwp...pppppp.r.reee.ee wwtpcccpppppppppceee httcc/ppprfffeesss.speeeeps files/ paper/
pdf
Yuheng Du, Qianhui Huang, Cedric Arisdakessian, and Lana X Garmire. Evaluation STAR and kallisto on single cell RNA-seq data aligument. G3 Genes Genomes| Genetics, 10(5):1775 1783, May 2020, doi:
URL https://doi.org::
1534//3. 120. 401160
Flavia Esposito. A review initialization methods for nonnegative matrix factorization: Towards omics
data experiments. Mathematics, 99), 2021. ISSN 2227-7390. doi: 1..3390/math90910000 URL https:
www..mdpi.cco//22777777077 wtt.://///2/2790/////////// //www.mdpi.c com
Vojtäch Franc, Václav Hlaváč, and Mirko Navara. Sequential coordinato-wise algorithm for the non-negative least squares problem. In Computer Analysis of Images and Patterns, pp. 407-414. Springer Berlin Hei- delberg. 2005. doi: 10.1007/11556121 50. URL
20

---

Published in Transactions on Machine Learning Research (04/2025)
Renaud Gaujoux and Cathal Seoighe. A flexible R package for nonnegative matrix factorization. BMC Bioinformatics, 11(1), July 2010. ISSN 1471-2105. doi: 10.1186/1471-2100--1-333 http:/t/://dd doi. org/10. 098/100.1186/17122205---111
Renaud Gaujoux and Cathal Seoighe. NMF: Algorithms and framework for nonnegative matrix factorization (NMF), 2023. URL https://CRAN, R- project ptt.a//o.........ooo.o..........
org/package-NFF
Gal Gilad, Itay Sason, and Roded Sharan. An automated approach for determining the number of com- ponents in non-negative matrix lactorization with application to mutational signature learning. Machine Learning: Science and Technology, 211)0015, December 2020. doi: 01188/2262-21511 2153/abc60a, URL
https://doi. org. wtt.aorrr//...//888.81 2153//22213/// /a//c0000
Nicolas Gillis. Sparse and unique nonnegative matrix factorization through data preprocessing. Journal of
Machine Learning Research, 13(108))3349 3386, 2012. URL http:://jml...rg/paapppaaaa wwt.elrrrrrrraapaeeeaeeeeeee wttp.///pppprrr///////ll/s12
html
Nicolas Gillis. Nonnegative matrix factorization. January 2020. doi: 10.1137/19978161197761
URL
Nicolas Gillis and François Glincur. Accelerated multiplicative updates and hicrarchical ALS algorithms for nonnegative matrix factorization. arXiv, 2011. doi: 00.48550/ARXIV.1107..0 URL https://arxix: org/ .194
T. R. Golub, D. K. Slonim, P. Tamayo, C. Huard, M. Gaasenbeek, J. P. Mcsirov, H. Coller, M. L. Loh, J. R. Downing, M. A. Caligiuri, C. D. Bloomfield, and E. S. Lander. Molecular classification of cancer: Class discovery and class prediction by gene expression monitoring. Science, 286(5439):531 537, October 1999. doi: 2..1266sccinnce 286.5339.666..5 URL httpps://oo. wtt.:/oo......//.c.ee.eoeoeo
286. .519..511
Lukas Grossberger, Francesco P. Battaglia, and Martin Vinck. Unsupervised clustering of temporal patterns in high-dimensionall neuronal ensembles using novel dissimilarity measure. PLOS Computational Bi- ology, e July 2018. doi: 10.1371/journal, 100.6.0288..... 1006283. URL https://doi. htttp:://ooi.or//1..55 http./..1//111111111111111
journal. pebraal.pbbi..m pcbi.cpii...... 1006283
D. Guillamet, B. Schiele, and J. Vitria. Analyzing non-negative matrix factorization for image classification. In 2002 International Conference on Pattern Recognition, volume 2, pp. 116-119 vol.2, 2002. doi: 10. 1109/TCCPR.2002.1445551
Ali Hassani, Amir Iranmanesh, and Najme Mansouri. Text mining using nonnegalive matrix factorization and latent semantic analysis, 2019. URL
04700
Cho-Jui Hsieh and Inderjit S. Dhillon. Fast coordinate descent methods with variable selection for non- nogative matrix factorization. In ACM SIGKDD International Conference on Knowledge Discovery and Dula Mining (KDD), aug 2011.
Lawrence Hubert and Phipps Arabie. Comparing partitions. Journal of Classification, 2(1):193 218, De- cember 1985. doi: 10. 10.11077/f000000 bf0199075.. URL https://ooaa.rr////.110///000000 1007/b/0990000
Yu Ito, Shin ichi Oeda, and Kenji Yamanishi. Rank selection for non- negatiee matrix factorization with normalized maximum likelihood coding. Proceedings of the 2016 SIAM International Conference on Dala Mining. Society for Industrial and Applied Mathematics, June 2016. doi: 10.1137/1.97816119778888.
81. URL httph:ttt//oooo otwpei..rg11... 1137/1 9T111//1.9771111744..88.
Bhargav Kanagal and Vikas Sindhwani. Rank selection in low-rank matrix approximations: : A study of cross-validation for nmfs. Proc Conf Adv Neural Inf Process, 2010.
Hyunsoo Kim and Hacsun Park. Sparse non-negative matrix factorizations via alternating non-ncgativity- constrained least squares for microarray data analysis. Bioinformatiis, 23(12):1495-1500) May 2007. doi: 10.1093/ 10.1093/bioinfoe /ttp///iforoaaiis//mmmaa3.. httpa///////////////////////////
21

---

Published in Transactions on Machine Learning Research (04/2025)
Jingu Kim and Haesun Park. Sparse nonnegative matrix factorization for clustering. Technical re- port, Georgia Institute of Technology, 2008. URT https://facultty:: ct..ccc...... www.eltttt....cccc......oa.e.eee wtwp.heeah//aaa/pparrrr..
GT CSE-08-01.pd
Keigo Kimura, Yuzuru Tanaka, and Mineichi Kudo. A fast hierarchical alternating least squares algorithm for orthogonal nonnegative matrix factorization. In Dinh Phung and Hang Li (eds.), Proceedings of the Sioth Asian Conference on Machine Learning, volume 39 of Proceedings of Machine Learning Research, pp. 129-141, Nha Trang City, Victnam, 26-28 Nov 2015. PMLR. URL httpss:/pproceednggss v39/kimural4.htr
press/
Vladimir Yu Kiselev, Tallulah S. Andrews, and Martin Hemberg. Challenges in unsupervised clustering of single-cell RNA-scq data. Nature Reviews Genetics, 20(5):273 282, Jamuary 2019. ISSN 1471-0064. doi:
10.1035/<41576-01--0088- URL http:://x..dd doi. wwwpb/ddoo.ooog....ooo http.//..........111///11101-000 110110011000//41777-------------
Kamran Kowsari, Donald F. Brown, Mojtaba Heidarysafa, Kiana Jafari Meimandi, Matthew S. Gerber, and Laura E. Barnes. HDLTex: Hierarchical deep learning for text classification. In 2017 16th IEEE International Conference on Machine Learning and Applications (ICMLA), pp. 364-371, 2017. doi: 10. 1109/ICMLA.2017.00
Kamran Kowsari, Donald E. Brown, Mojtaba Heidarysafa, Kiana Jafari Meimandi, Matthew Gerber, and Laura E. Barnes. Web of science dataset, March 2019. URL https://dataa.mendollo. www.eenooooll..ommmtttseeeeeeeee com/datasets/
9rw3vkcfy4/6
Ken Lang. 20 newsgroups, 1997. URL htte://ggonn.ccmm/aaaaa///////gp
Amy N. Langville, Carl D. Meyer, Russell Albright, James Cox, and David Duling. Algorithms, initializations, and convergence for the nonnegative matrix factorization, 2014. URL bttps://arxiv.org abs/ 1407. 7299
Daniel Lee and H. Sebastian Soung. Algorithms for non-negative matrix factorization. In T. Leen, T. Dietterich, and V. Tresp (eds.), Advances in Newral Information Processing Systems, vol- ume 13. MIT Press, 2000. URL
f9d1152547c0bdee116 fa0111244700ae133300775560022622
pdf
Daniel D. Lee and H. Sebastian Seung. Learning the parts of objects by non-negative matrix factorization. Nature, 401 (6755)))788777) -791, October 1999. doi: 10.1038/44565.
Chih-Jen Lin. Projected gradient methods for nonnegative matrix factorization. Newral Computation, 19 (10):2756 2779, October 2007. doi: n..1622/eeo.2007.111 19.12..5. URL
neco 200...007..00 19. 19.10.27.6..
Xihui Lin and Paul C. Boutros. Optimization and expansion non-negative matrix factorization. BMC Bioinformaties, 21(1), January 2020. doi: (0.1186/912859-019 812859-019-3312-1
URL httpps////// http:///doo.ogg//...1111 org/10.1186/
Haifeng Liu, Zhaohui Wu, Xuelong Li, Deng Cai, and Thomas S. Huang. Constrained nonnegative matrix factorization for image representation. IEFE Transactions on Pattern Analysis and Machine Intelligence, (4(7))1299-1331) 2012. doi: 10.100//PAAM12201122..
Laura Muzzarelli, Susanne Weis, Simon B. Eickhoff, and Kaustubh R. Patil. Rank selection in non-negative matrix factorization: systematic comparison and new MAD metric. In 2019 International Joint Con- ference on Neural Networks (IJCNN), pp. 1 1-8, 2019. doi: 10.1109/1JCNN..C 2011
Art B. Owen and Patrick O. Perry. vi-croos-vllidation of the SVD and the nonnogative matrix factorization. The Annals of Applied Statistics, 3(2), June 2009. doi: 10.1214 111121/4--a0 101/088aosss aoas227. URL https: http:///ooooo////.//////// 1214/08- aoas227
22

---

Published in Transactions on Machine Learning Research (04/2025)
V. Paul Pauca, Farial Shahnaz, Michael W. Berry, and Robert J. Plemmons. Text mining using non- negative matrix factorizations. In Proceedings of the 2004 SIAM International Conference on Data Mining, Society for Industrial and Applied Mathematics, April 2004. doi: 10.1137/1978111117774444 URL https:///oi.ors www.ooo..ogo....... 10. 1137/1 978111197771111 45
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapcau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825 2830, 2011.
Rafael Rosales, Rodrigo D Drummond, Renan Valieris, Emmanuel Dias-Neto, and Israel T da Silva. signeR: an empirical Bayesian approach to mutational signature discovery. Bioinformatics, 33(1): 16, September 2016. ISSN 1367-4811. doi: Lt..B///boiffoomatcccss///////// btt. URL http://dxx doi oww..///oo...o....//1
1093/biolnffrrmatccc/bbbaa
Mikkel N. Schmidt, Ole Winther, and Lars Kai Hansen. Bayesian non-negative matrix factorization. In Independent Component Analysis and Signal Separation, pp. 540-547. Springer Berlin Heidelberg, 2009. doi: 1110007977-3-612--00
68. URL hattpa://don org/ 10.
Gregory W. Schwartz, Yeqiao Zhou, Jelena Petrovic, Maria Fasolino, Lanwei Xu, Sydney M. Shaffer, Warren S. Pear, Golnaz Vahedi, and Robert B. Faryabi. Toomanyeells identifies and visualizes re- lationships of single-cell clades. Nature Methods, 17(4):405-44) 413, March 2020. ISSN 1548-7105. doi:
URL http://dx doi wwt.:/ooooo.oog/.//// 10. 2tt////42122222-00070
Aristeidis Sotiras, Jon B. Toledo, Raquel E. Gur, Ruben C. Gur, Theodore D. Satterthwaite, and Christos Davatzikos. Patterns of coordinated cortical remodeling during adolescence and their associations with functional specialization and evolutionary expansion. Proceedings of the National Academy of Sciences, 114(13) 3527-3333333333 1620928114
March 2017. doi: 10. (0.1073/pnas 162092862092114.. URL
V. Y. F. Tan and C. Févotte. Automatic relevance determination in nonnegative matrix factorization. In Proc. Workshop on Signal Processing with Adaptative Sparse Structured Representations (SPARS), St- Malo, France, Apr. 2009. URL 10. 1109/TPAMM 2012. 240.
Bosiljka Tasic, Vilas Menon, Thuc Nghi Nguyen, Tae Kyung Kim, Tim Jarsky, Zizhen Yao, Boaz Levi, Lucas T Gray, Staci A Sorensen, Tim Dolbeare, Darren Bertagnolli, Jeff Goldy, Nadiya Shapovalova, Sheana Parry, Changkyu Lee, Kimberly Smith, Amy Bernard, Linda Madisen, Susan M Sunkin, Michael Hawrylycz, Christof Koch, and Hongkui Zeng. Adult mouse cortical cell taxonomy revealed by single cell transcriptomics. Nature Neuroscience, 19(2):335 346, January 2016. ISSN 1546-1726. doi: 10.1038 nn. 4216. URL http:://sss doi. org/ 1038/ nn. 4216
Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Ev- geni Burovski, Pearu Peterson, Warren Weckessor, Jonathan Bright, Stéfan van der Walt, Matthew Breut, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, îhaa Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261 272, 2020. doi:
Yu-Xiong Wang and Yu-Jin Zhang. Nonnegative matrix factorization: A comprehensive review. TEEE Transactions on Knowledge and Data Engineering, 25(6):1336 1353, 2013. doi: 10.1109/TKDE:22222.51.
Svante Wold. Cross-validatorry estimation of the number of components in factor and principal components models. Technometrics, 20(4):397 405, November 1978. doi: 10.1080/00401100.11973 https:///oi. org/ 10. 11tp///...////0000000111 1978. 111.778.1189999 10489693
URL
Zhijin Wu and Hao Wu. Accounting for cell type hierarchy in evaluating single cell RNA-seq clustering. Genome Biology, 21(1), May 2020. ISSN 1474-760X. doi: 10.1186/s13059-020 x URL http:
7/dx. wwwp:.doo.org/1..1 1186/ (1118//1100 -13059- -200000900000002-2030000
23

---

Published in Transactions on Machine Learning Research (04/2025)
Sayaka Yamauchi, Masanori Kawakita, and Jun'ichi Takeuchi. Botnet detection based on non-negative matrix factorization and the MDL principle. In Newral Information. Processing, pp. 400 100-409.00 Springer Berlin Heidelberg, 2012. doi: 10.10077977--3-6223455000: 48. URL https://doi.co /ttp///....rg11.100077//
378-3-642- 5 48
Xiaohui Yang, Wenming Wu, Xin Xin, Limin Su, and Liugen Xuc. Adaptive factorization rank selection- based NMF and its application in tumor recognition. International Journal of Machine Learning and Cybernetics, 12(9))2673 2691, May 2021. ISSN 1088-00X. doi: 10.1007/9133422-211-0133 1. URL http: doi. (rdd..o.....//0001000///02020000
Marinka Zitnik and Blaz Zupan. Nimfa: A Python library for nonnegative matrix factorization. Journal of Machine Learning Research, 13:849 853, 2012.
24


Reviews:
[{"review_id": "pfwD1k6DaO", "reviewer": "Reviewer_XcKY", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n1. The paper is easy to read and follow.\n2. The paper studies an interesting problem of the matrix factorization community.\n3. The proposed approach is novel.\n\n\n### Weaknesses\n\n1. **The experiments are based on the concept of a \"true rank\" which I consider problematic.** The authors state in Section 4 that they define the \"true rank\" of a dataset as the number of classes of the underlying classification dataset and then base their evaluation on that \"true rank\". I think this is highly problematic. NMF embeds data within a cone and the rank $k$ can be seen as the number of vertices. While it would be nice that every class has its own vertex (corner), this is not necessarily true. Consider handwritten digits. One could hope that 10 is the perfect rank as there are 10 digits, but there are different variants to write digits, thus, one could also argue for a $k$ that is higher than 10. Another example would be where a class is not a single mode but rather multi-modal.\n2. **The experiments are not convincing.** Even if we believe that there is a \"true rank\" which is the number of classes of a classification dataset, I do not consider the experiments convincing. Consider any table. The proposed method MCI-RSIC is not really better than any other baseline. While it sounds nice that the proposed method can suggest multiple ranks, which one should be chosen? I do not often see the \"true rank\" within the suggestion. Besides, one could also argue that the elbow method is able to suggest multiple ranks.\n3. **Time is crucial factor which was not at all evaluated.** Assume a practitioner is interested in a new method for rank suggestion. A critical aspect apart from the performance is time. The proposed method MCI-RSIC relies on many different initializations per potential rank. This is time-intensive. Thus, computation time should be reported. It could be that the suggested method is factors slower than all baselines which often also give reasonable rank estimates. Is the effort worth it? The reader does not know.\n4. **The initialization of NMF is very important and it is insufficiently discussed and evaluated.** As mentioned before, many different initializations per potential rank. It is known that the performance of NMF (like many other methods) heavily relies on the initialization. I am missing a discussion of the state of the art when it comes to initialization a NMF model. Would a better initialization or a more diverse set of initializations help? It remains unclear.\n5. **An ablation study is missing.** Continuing the previous point regarding the initializations (would initializations other than uniform be better?), it would be interesting to see whether 100 initializations are really necessary. It is unclear how much the rank suggestion relies on the number of initializations.\n6. **The concept of \"islands of stability\" is weird.** In Figure 2, 29 (and 38) is no \"island of stability\" although 19 and 23 are. Why? Because the bump is not as big? But then in Figure 3, 21 is an \"island of stability\".\n7. **The order of presentation is not always ideal.** Experiment-specific choices should reside in the experiment section but are presented throughout the paper, see Section 2.1 (how to optimize NMF), Section 2.3.3 (evaluation criteria), Section 2.3.4, and Section 3.\n\n\n### Questions\n1. What are the rows and columns in Figure 1? What is the range of the axes (0-1200)? If just the first 400 images are meant, please state 400 in the caption and drop the axes ranges.\n2. I have a problem with the statement \"We define the 'true rank' of a dataset as the number of classes in the underlying dataset.\" Please see my weakness above. One can easily define new datasets where the number of dimensions, classes, and optimal rank is arbitrary, meaning the number of classes is different from the optimal rank. Why should the number of classes be a good indicator for the \"true rank\"?\n3. If you force the optimization to run 100 iterations, why do you also set a tolerance? Either you always ran 100 iterations or you stopped when the tolerance criteria is met, but not both. Can you clarify?\n4. How can I interpret the seed 123456789? Was it set once before the 100 initializations? You could have also set a seed of $r$ for the $r$-th initialization.\n5. In Section 6, RAM constraints were mentioned. I wonder how much RAM there was. Runtimes of 30 days seem unrealistically long for such small problems. Is a single run meant, meaning one optimization for one initialization?\n6. Figure 2: If ranks 1-5 all have similar MCI values, why isn't rank 1 preferred over rank 5? Isn't the simpler model always better?\n7. Considering the result tables, e.g., Table 1, why should the reader choose the proposed approach given that the true rank is 2, all baselines suggest 2-5 and MCI greater than 5? The elbow seems to do the best job here.\n8. In Figure 2, why are 29 and 38 not considered as islands of stability?\n\n\n### Minor Comments\n- Abstract: Typo in \"Intial\"\n- Abstract: Why is Mean Coordinatewise Interquartile Range abbreviated with MCI and not MCIR?\n- Abstract: Why is there \"existing\" in \"current state-of-the-art existing rank\", if it is the current SOTA, it has exist.\n- Introduction, first paragraph: The phrase \"complex datasets\" can be misleading (datasets consisting of complex numbers).\n- Introduction, first paragraph: What are examples for \"considerable limitations\"?\n- Page 2, first paragraph: The references for \"Other methods\" are missing.\n- Page 2, first paragraph (occurs more often): cross-validation-based\n- Page 2, first paragraph: It sounds like the preference for lower ranks is a bad thing or not intentional, but in doubt we always prefer the simpler thing. Think about Occam's razor.\n- Page 2, second paragraph (occurs more often): Please remove \"given\" and capitalize section when referring to a specific one.\n- Section 2.1: \"Although a variety methods\" is missing a word.\n- Section 2.1: Arbitrary words are capitalized. Why?\n- Section 2.1: Space before full-stop: \"in this study . Additionally\"\n- Section 2.1: Two objectives are mentioned but then you choose a third. This is weird. The Frobenius norm is not a single Euclidean distance as suggested in the text.\n- Section 2.1: SCD updates: the equations only work after initialization\n- Section 2.1: The clustering aspect of NMF drops out of no-where. This should be rephrased.\n- Section 2.2: \"There are typically numerous equally as good solutions (e.g. factorizations for which the Frobenius norm of the residual are equally small)\". Remove \"as\", use \"norms\" and \"residuals\", the norms could also be equally \"big\"!\n- Section 2.2: Full-stop missing: \"2008) Third,\"\n- Section 2.3.1: \"unless choosing to run on a limited range of ranks\" The range of ranks is always limited!\n- Section 2.3.2: \"randomly in half\"\n- Section 2.3.2: Why is Rand capitalized (while adjusted and index are not). Same with AIC in Section 2.3.3.\n- Section 2.3.3: Wrong citation style for Yamauchi et al. (2012), and Tan & Fevotte (2009).\n- Section 2.3.4: References are missing (MADImput, etc).\n- Sections 2.3.4, 4.1, 5.3, 5.5: Many wrong citation styles.\n- Section 2.4: The meaning of $\\sim$ should be introduced earlier (at least one sentence earlier).\n- Section 2.4: Equation (3), please use \\eqref\n- Page 7, last equation: Please use \\text{opt} and \\text{init} as before. If opt refers to optimized, I suggest to use NMF, since opt could be interpreted as optimal, which is not true!\n- Section 3.2: Please use \\cdot for $\\text{vec}(\\cdot)$.\n- Section 3.2: Rephrase \"may now\" \n- Section 4: If a dataset is taken from scikit-learn, please name scikit-learn and cite the paper once. The reader might be familiar with scikit-learn but less so with the last name of the first author of this publication.\n- Section 5: There are a lot of tense mismatches. Either, experiments happen as the reader reads the manuscript, or they all happened in the past, but the tenses are mixed.\n- Page 10, first paragraph: The Frobenius norm is not defined in Equation (1)!\n- Page 10, first paragraph: How much memory and cores does the other machine have? The information is insufficient.\n- Section 5.3: \"the the\"\n- Section 5.3: State numpy. However, this sentence (package for shuffling) can also be deleted as it has no useful information.\n- Section 5.3: Why is now the median computed when before, the average was computed (Section 5.1)?\n- Section 6: A result table is mentioned but it is unclear which one. So far, no table appeared.\n- Section 6: I do not understand why Section 5 states that $k_{\\text{min}}=2$ but in Section 6 it changes to $k_{\\text{min}}=1$ and the figures start from zero. The figures should start from $k_{\\text{min}}$ and end with $k_{\\text{max}}$. Axis labels should be provided (not only in the caption) and a light grid would make it easier to read.\n- Figure 6 is way too big.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Non-negative matrix factorization decomposes a matrix with non-negative elements into a factor and weight matrix, both of which are also non-negative. Like many other unsupervised learning methods such as $k$-means clustering, NMF relies on the specification of a rank $k$, i.e, the latent dimension of the factorization. One well-known example for choosing the rank is the elbow method, where the objective function value is plotted against the rank. The suggested rank is the first rank after which the improvement in terms of objective function value is strongly diminishing.\n\nThis paper proposes a new method called Mean Coordinatewise Interquartile Range Residual Sensitivity to Initial Conditions (MCI-RSIC) which aims at suggesting multiple suitable ranks instead of a single rank. The idea is to compute, per considered rank $k$, a number of correlated initializations, compute NMF per initialization, compute the residuals per NMF (per initialization), and finalluy to compute the mean coordinatewise interquartile range. After plotting the proposed metric against the potential ranks $k$, so-called \"islands of stability\" can be identified where the metric changes.\n\nThe authors compare their proposed approach on a number of (classification) datasets against baseline methods for rank suggestion for NMF.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please respond to my questions and weakness points and implement the changes.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "I8rQAi1FIT", "reviewer": "Reviewer_sp8L", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\nThe description on previous work on Rank Determination is appreciated, and allows the reader to understand the current research landscape on Rank Determination, as well as the strengths and weaknesses of each method. Additionally, the diverse set of datasets in the numerical analysis is also useful and interesting.\n\n**Weaknesses**\n\nThis work seems to me as an interesting proof of concept for a novel Rank Determination method. However, I feel that the main weakness is the completeness of the paper itself. I am unsure if the scope of study is enough for publication as a journal paper. I have written some comments below and hope they are useful.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work proposes a new method for determining ranks in Non-negative Matrix Factorization problems. The main contribution is the proposal of a novel approach (RSIC) that determines a selection of ranks as opposed to one rank. The numerical experiments on a large selection of data suggest that the proposed method's performance is comparable to existing methods.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1) **Computational Complexity:** This work emphasizes the complexity burden that methods such as Cross Validation present. It also presents RSIC as a \"scalable\" and \"generalizable\" method without any study of the complexity of the method. Even a numerical analysis of the complexity (such as comparing computation times with other methods) would be appreciated.\n\n2) **Usefulness of Approach:** It would be beneficial for the author(s) to further discuss why knowing a range of ranks may be useful, or introducing (even briefly) a specific application where one may be more interested in knowing a collection of ranks of interest, or where knowledge of the \"island of stability\" can be leveraged.\n\n3) **Choosing the Rank:** The question of which rank is ultimately chosen from a collection of computed ranks is, in my view, an important open question and a crucial part of this work. In Table 2, MCI-RSIC performs comparably to Elbow, even when the latter \"undershoots\". MCI-RSIC provides a collection of ranks but the wrong choice of rank will also lead to undershooting. I appreciate that the author(s) recognized this ad future work but some discussion or insights must also be included in this work.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "9r3VpK2CMT", "reviewer": "Reviewer_iw4J", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n\n1. The paper investigates an interesting problem, that is exploring the relationship between the stability of non-negative matrix factorization and the target rank.\n2. The paper proposes a new metric, MCI, that could potentially capture multiple ranks with relatively good statistical stability in NMF.\n3. The paper conducts a large set of numerical simulations to demonstrate how to get the potential \"good\" ranks in NMF with MCI metric.\n\nWeakness:\n\n1. The proposed solution, identify \"island of stability\" using the MCI metric, is not rigorously defined. To my best knowledge, the criteria of finding the \"island of stability\" is only verbally explained without any rigorous mathematical definition in Sec 3. Moreover, the verbal explanation is \"island of stability\" exhibits \"diminished variance\" or \"minimal compared to nearby ranks\", however, there is no clear definition of how much decrease is diminishing, and how many neighbour ranks count as nearby. Later on, the criteria seem to change to the point before large spread increase. The criteria need to be carefully defined and discussed.\n\n2. To identify \"good\" ranks, the authors are finding ranks that follows the vaguely defined definition. However, note that the MCI metric is an empirical average of residuals. Therefore, error bars should be provided in Figure 2-5 to convince readers that the identified ranks are not simply due to numerical errors.\n\n3. The authors did not provide any explanation behind the \"goodness\" of ranks besides the simple argument of the more stable the better. This is partly reflected by that fact that the \"good\" ranks identified by MCI metric are usually not the same as the \"true rank\". \n\nOverall: The paper, at this moment, serves more as a vaguely-written manual for a new approach in identifying stable ranks in NMF, rather than a scientific paper that provides empirical and theoretical justification for the proposed approach. Although the quantity of numerical simulations is large, the quality of the simulations needs improvement. A smaller set of simulations that examines different aspects of the algorithm with more nuanced discussion, in my opinion, is much better than the current simulation section.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper argues that non-negative matrix factorization is unstable due to its random initialization + optimization nature. Moreover, the authors argue that there could exist several \"good\" ranks for a given dataset, and the \"goodness\" is reflected by the spread of solution residuals under a certain rank with random initializations. To this end, the authors proposed to find the \"good\" ranks with MCI metric over multiple random initializations, where lower MCI metric indicates the optimal NMF solutions are less unstable. The authors provides several simulations to demonstrate how to find the good ranks for different datasets.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please see weakness.\n\nMinor: \n\n1. Figure 1 is not so informative as it just seems errors are scattered throughout the datasets, which is a common phenomenon. The authors could explain it more in detail, preferably with some mathematical language.\n\n2. There are some small typos throughout the paper, e.g., in Sec 2.2, though it (\"is\" missing) largely domain-dependent, the footnotes 3 in Table 1 points to nowhere.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers and authors engaged in a robust and fruitful discussion. The comments by the reviewers and changes to the manuscript resulted in a much improved submission. There was substantial discussion around the use of evaluative words like "best" and comparison to a "true" rank. These concerns are well founded and the authors made efforts to adjust the claims to rest firmly on the foundation of the evidence. However, the community norms also impinge on the word selection - the authors have little control over those expectations in the short run, though one expects that the field will converge to appropriate descriptions in the long run. Within the bounds of the expectation of the journal the authors have responded to the most important comments. Of course, more work may be done to test the methodology in different ways. The paper seems to meet the standards of the evaluation criteria.
Decision: Accept (poster)