A New Family of Poisson Non-negative Matrix Factorization Methods Using the Shifted Log Link

Eric Weine; Peter Carbonetto; Rafael A. Irizarry; Matthew Stephens

cs.LGJan 9, 2026

A New Family of Poisson Non-negative Matrix Factorization Methods Using the Shifted Log Link

Authors: Eric Weine, Peter Carbonetto, Rafael A. Irizarry, Matthew Stephens

Organizations: Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA, USA · Department of Human Genetics, University of Chicago, Chicago, IL, USA · Department of Data Science, Dana-Farber Cancer Institute, Boston, MA, USA · Department of Statistics and Department of Human Genetics, University of Chicago, Chicago, IL, USA

Abstract

Poisson non-negative matrix factorization (NMF) is a widely used method to find interpretable "parts-based" decompositions of count data. While many variants of Poisson NMF exist, existing methods assume that the "parts" in the decomposition combine additively. This assumption may be natural in some settings, but not in others. Here we introduce Poisson NMF with the shifted-log link function to relax this assumption. The shifted-log link function has a single tuning parameter, and as this parameter varies the model changes from assuming that parts combine additively (i.e., standard Poisson NMF) to assuming that parts combine more multiplicatively. We provide an algorithm to fit this model by maximum likelihood, and also an approximation that substantially reduces computation time for large, sparse datasets (computations scale with the number of non-zero entries in the data matrix). We illustrate these new methods on a variety of real datasets. Our examples show how the choice of link function in Poisson NMF can substantively impact the results, and how in some settings the use of a shifted-log link function may improve interpretability compared with the standard, additive link.

Figures & tables

Figure 1: Plots of the link function $g(\lambda;c)=\alpha_{c}\times\log(1+\lambda/c)$ for various values of $c$ .

Model	Log-likelihood Complexity
Standard Poisson NMF	$\mathcal{O}((\omega+n+p)K)$ ( Carbonetto et al., 2021 )
Poisson GLM-PCA	$\mathcal{O}(npK)$ ( Weine et al., 2024 )
Log1p Poisson NMF Exact	$\mathcal{O}(npK)$
Log1p Poisson NMF Approximate	$\mathcal{O}\!\left((\omega+n+p)K+(n+p)K^{2}\right)$

Table 1: Computational complexity of log-likelihood for Poisson matrix factorization variants on sparse data matrices (size

n\times p

with

\omega

non-zero entries).

Figure 2: Scaling of computational complexity of equations ( 12 ) and ( 15 ) with data characteristics and choice of $K$ . The y-axes are the ratio of the computational complexity of the specified calculation relative to the complexity of computing the log-likelihood of the standard Poisson NMF model, where $1$ indicates the two calculations have the same complexity. (A) Scaling with respect to $n$ with all other variables fixed (note that this is equivalent to scaling with $p$ ). (B) Scaling with $K$ with all other variables fixed. (C) Scaling with sparsity $\left(\left[1-\frac{\omega}{np}\right]\cdot 100\%\right)$ with all other variables fixed.

Figure 3: Per-entry likelihood ratio (the likelihood of the data when optimized with the approximate objective divided by the likelihood when optimized with the exact objective, raised to the power $1/np$ ) of factor models fit with $K=5$ and varying settings of $c$ . Points show means over simulation replicates ( $30$ per generative regime; $100$ for $c=10^{-3}$ ) and error bars show $95\%$ confidence intervals. All data are generated with $n=p=500$ and $K=5$ . (A) Likelihood ratios when data are generated from the log1p model with $c=10^{-3}$ . (B) Likelihood ratios when data are generated from the log1p model with $c=1$ . (C) Likelihood ratios when data are generated from the log1p model with $c=\infty$ .

Figure 4: Toy example results. (A) True feature rates in each of the four sample groups. (B) True sample scores, and estimated sample scores ( $\mathbf{L}$ ) for log1p NMF with $c=\infty$ and $c=1$ . (C) Distributions of estimated factor values (columns of $\mathbf{F}$ ) for the two change factors, with features classified by their baseline level and by whether the change factor increases them.

Figure 5: Results of mis-specification simulation: Thin lines show individual simulation replicates, thick lines their mean, and the open circle marks $c_{\textrm{fit}}=c_{\textrm{true}}$ . (A, B) Mean Hoyer sparsity of the columns of $\mathbf{L}$ and of $\mathbf{F}$ ; the dashed line shows the same quantity for the true loadings and factors. (C) Peak-to-mean ratio $\operatorname{pm}(\mathbf{B})$ of the fitted model; the dashed line shows the value predicted from the true rates alone — $\operatorname{pm}(\log(1+\boldsymbol{\Lambda}/c))$ at finite fitting $c$ , and $\operatorname{pm}(\boldsymbol{\Lambda})$ at $c=\infty$ — averaged over replicates. (D) Per-entry likelihood ratio relative to the same replicate’s fit at $c_{\textrm{true}}$ , i.e. $\exp\{[\ell_{\textrm{log1p}}(c)-\ell_{\textrm{log1p}}(c_{\textrm{true}})]/np\}$ ; the dashed line marks $1$ .

Figure 6: Combined figure and table for the MCF-7 analysis. (A-B) Visual representation of estimated $\mathbf{L}$ matrices for log1p NMF with $c=1$ and $c=\infty$ . Each column represents a row of $\mathbf{L}$ , where each color corresponds to a column of $\mathbf{L}$ . (C-D) Scatterplots of factors $2$ and $3$ . Each point corresponds to a single gene (row of $\mathbf{F}$ ). Points are colored based on results of differential expression using DESeq2 ( Love et al., 2014 ) . Points in green have Benjamini-Hochberg ( Benjamini and Hochberg, 1995 ) adjusted p-values $<0.01$ and log2FC $>1$ in both the RA and TGF- $\beta$ groups. Points in orange meet these conditions in only the RA group, points in blue meet these conditions only in the TGF- $\beta$ group, and points in grey meet these conditions in neither group. (Table): Top $10$ genes in each column of $\mathbf{F}$ from the fitted models. Below each gene list, in upright type, is the most significantly enriched Gene Ontology ( Ashburner et al., 2000 ; The Gene Ontology Consortium, 2023 ) biological process term among the factor’s top $30$ genes, computed by a hypergeometric test against all genes in the dataset ( Boyle et al., 2004 ) , with Benjamini–Hochberg adjusted $p<0.05$ .

Figure 7: Combined figure and table for the cell type associated factors of the Pancreas analysis. (A-B) Visual representation of estimated $\mathbf{L}$ matrices for the log1p NMF with $c=1$ and $c=\infty$ , grouped by cell type. Each column represents a row of $\mathbf{L}$ , where each color corresponds to a column of $\mathbf{L}$ . (Table) Top genes for the factors of each model, excluding treatment associated factors. GO terms are computed as described in the caption of Figure 6 .

Figure 8: Comparison between delta and gamma cells in single cell pancreas data. (A-B) Gene scores ( $\mathbf{F}$ ) from $c=1$ and $c=\infty$ models corresponding to the factor most associated with delta cells (k3) and the factor most associated with gamma cells (k13). Before plotting, each matrix $\mathbf{L}$ was normalized so that the maximum value of each column was $1$ , and $\mathbf{F}$ was scaled accordingly. (C) Mean expression of genes in delta and gamma cells.

Figure 9: Visualization of the treatment associated factors of the Pancreas analysis. (A-B) Cell scores grouped by treatment. (C-D) Cell scores grouped by cell type. (Table) Top genes for the factors of each model, excluding cell type associated factors. GO terms are computed as described in the caption of Figure 6 .

Figure 10: Combined figure and table for the BBC analysis. (A-B) Visual representation of estimated $\mathbf{L}$ matrices for log1p NMF with $c=0.001$ and $c=\infty$ , grouped by document type. Each column represents a row of $\mathbf{L}$ , where each color corresponds to a column of $\mathbf{L}$ . (Table) Top words for each factor.

Figure 11: Summary of model fits across values of $c$ for real data applications. Spearman correlations are calculated pairwise between all factors (columns of $\mathbf{F}$ ) and averaged. Sparsity values are calculated for each loading / factor (columns of $\mathbf{L}$ or columns of $\mathbf{F}$ ) and then averaged. (A) Correlation and sparsity values for the MCF-7 dataset with $K=3$ . (B) Correlation and sparsity values for the pancreas dataset with $K=13$ . (C) Correlation and sparsity values for the BBC dataset with $K=10$ .

Appendix figures & tables9 assets

Supplementary material from the paper’s appendix.

Appendix

Figure F.1: Progress of the exact and approximate (Chebyshev) optimization approaches on three real datasets, fitting log1p NMF with $c=1$ from the rank-1 initialization. Each panel shows, as a function of wall-clock time (including initialization time), the distance between the exact log-likelihood of the current iterate and the best exact log-likelihood achieved by either approach on that dataset; a curve reaching $0$ attained the best solution found. Distances are plotted on a $\log(1+x)$ scale, and the time axes of panels (B) and (C) are on the same scale, so that $0$ (the initialization) remains representable. (A) MOCA. (B) Pancreas. (C) BBC.

Figure G.1: Celltype associated factors for the log1p NMF model fit by maximizing the exact log-likelihood (A), as well as with the Chebyshev approximation approach (B), and by fitting Frobenius NMF to the log1p transformed counts (C).

Figure G.2: Treatment associated factors for the log1p NMF model fit by maximizing the exact log-likelihood (A), as well as with the Chebyshev approximation approach (B), and by fitting Frobenius NMF to the log1p transformed counts (C).

Figure K.1: Number of iterations until convergence for the fits of the approximation accuracy experiment of Section 3.3 . For each generative regime and each finite value of $c$ used for fitting, points show the geometric mean of the number of iterations until convergence over simulation replicates ( $30$ per regime; $100$ for data simulated with $c=10^{-3}$ ), with error bars showing $95\%$ confidence intervals computed on the log scale.

Figure K.2: Stability of the $K=10$ BBC fits to initialization. For each $c\in\{0.001,1,\infty\}$ , the model was fit to convergence (optimizing the exact log-likelihood) from $30$ random initializations and from the rank-1 initialization used for the results in the main text. The random initializations are matched across values of $c$ : for each replicate we drew matrices $\mathbf{L}_{\textrm{raw}}$ and $\mathbf{F}_{\textrm{raw}}$ with i.i.d. $\textrm{Uniform}(0,1)$ entries, depending only on the replicate so that all values of $c$ share the same draws, and initialized each model at $\mathbf{L}_{0}=\sqrt{\tau_{c}}\,\mathbf{L}_{\textrm{raw}}$ , $\mathbf{F}_{0}=\sqrt{\tau_{c}}\,\mathbf{F}_{\textrm{raw}}$ , where $\tau_{c}>0$ was found numerically such that the Poisson rates implied by $\mathbf{B}_{0}=\tau_{c}\mathbf{L}_{\textrm{raw}}\mathbf{F}_{\textrm{raw}}^{\top}$ under that model’s link (i.e., $s_{i}g^{-1}(b_{ij};c)$ , and $\tau_{c}\mathbf{L}_{\textrm{raw}}\mathbf{F}_{\textrm{raw}}^{\top}$ directly for $c=\infty$ ) have mean equal to the mean of the observed counts. Light points show the exact log-likelihood of each converged random-initialization fit, centered at the mean over the $30$ random initializations for that value of $c$ ; error bars show $\pm 1$ standard deviation around this mean. (Table): mean and standard deviation of the exact log-likelihoods of the random-initialization fits for each value of $c$ .

Figure K.3: Additional results for Simulation 2. Thin lines show individual simulation replicates, thick lines their mean, and the open circle marks $c_{\textrm{fit}}=c_{\textrm{true}}$ . (A) RMSE of the fitted rates against the true rates, relative to the same replicate’s fit at $c_{\textrm{fit}}=c_{\textrm{true}}$ (log scale); the dashed line marks $1$ . (B) Mean absolute pairwise Spearman correlation between the columns of the estimated $\mathbf{F}$ (as in Figure 11 ); the dashed line shows the same quantity for the true factors.

Figure K.4: Structure plots for the MCF-7 dataset across various values of $c$ .

Figure K.5: Structure plots for the pancreas dataset across various values of $c$ .

Figure K.6: Structure plots for the BBC dataset across various values of $c$ .

Non-Negative Matrix Factorization

Non-negative matrix factorization (NMF) decomposes nonnegative data into additive, interpretable latent components for analysis and discovery.

See all 33 papers

MM-algorithms for traditional and convex NMF with Tweedie and Negative Binomial cost functions and empirical evaluation
Mar 10, 2026Elisabeth Sommer James, Asger Hobolth, Marta PelizzolaNon-Negative Matrix Factorization
Department of Mathematics, Aarhus University, Denmark
An Efficient Newton Algorithm for Nonnegative Matrix Factorization with the Kullback-Leibler Divergence
Jul 15, 2026Damien Lesens, Jérémy E. Cohen, Bora UçarNon-Negative Matrix Factorization Low-Rank Matrix Decomposition
ENS de Lyon · Univ. Lyon, INSA-Lyon, UCBL, UJM, CNRS, Inserm, CREATIS UMR 5220, U1294, F-69621, Lyon, France. · CNRS and LIP (UMR5668, Université de Lyon - ENS de Lyon - UCBL - CNRS - Inria), 46 allée d'Italie, Lyon, France.
Nonnegative Matrix Factorization in the Component-Wise L1 Norm for Sparse Data
Mar 31, 2026Giovanni Seraghiti, Kévin Dubrulle, Arnaud Vandaele +1Non-Negative Matrix Factorization Matrix Factorization
University of Mons, Rue de Houdain 9, 7000 Mons, Belgium, and Dipartimento di Ingegneria Industriale, Università degli Studi di Firenze, Viale Morgagni 40/44, 50134 Firenze, Italia. Member of the INdAM Research Group GNCS. · University of Mons, Rue de Houdain 9, 7000 Mons, Belgium.

Mar 10, 2026J/K move · Enter open · S save

MM-algorithms for traditional and convex NMF with Tweedie and Negative Binomial cost functions and empirical evaluation

Elisabeth Sommer James, Asger Hobolth, Marta Pelizzola

Department of Mathematics, Aarhus University, Denmark

Non-Negative Matrix Factorization

Open paper AskPDF arXiv

Non-negative matrix factorisation (NMF) is a widely used tool for unsupervised learning and feature extraction, with applications ranging from genomics to text analysis and signal processing. Standard formulations of NMF are typically derived under Gaussian or Poisson noise assumptions, which may be inadequate for data exhibiting overdispersion or other complex mean-variance relationships. In this paper, we develop a unified framework for both traditional and convex NMF under a broad class of distributional assumptions, including Negative Binomial and Tweedie models, where the connection between the Tweedie and the

β

-divergence is also highlighted. Using a Majorize-Minimisation approach, we derive multiplicative update rules for all considered models, and novel updates for convex NMF with Poisson and Negative Binomial cost functions. We provide a unified implementation of all considered models, including the first implementations of several convex NMF models. Empirical evaluations on mutational and word count data demonstrate that the choice of noise model critically affects model fit and feature recovery, and that convex NMF can provide an efficient and robust alternative to traditional NMF in scenarios where the number of classes is large. The code for our proposed updates is available in the R package nmfgenr and can be found at https://github.com/MartaPelizzola/nmfgenr.

Abstract

Non-Negative Matrix Factorization

Explore similar work

MM-algorithms for traditional and convex NMF with Tweedie and Negative Binomial cost functions and empirical evaluation

An Efficient Newton Algorithm for Nonnegative Matrix Factorization with the Kullback-Leibler Divergence

Nonnegative Matrix Factorization in the Component-Wise L1 Norm for Sparse Data