This work argues that the KL-NMF method has reached its limits and proposes to use instead the second-order Taylor expansion of the loss, leading to a Newton-type method which provably converges and competes favorably with state-of-the-art algorithms on a large variety of datasets.
Abstract
Nonnegative Matrix Factorization (NMF) is a fundamental tool in unsupervised learning, which approximates a nonnegative matrix by the product of two low-rank nonnegative factors. The Kullback-Leibler (KL) divergence is best suited to measure the data to model discrepancy when the decomposed data sample follows a Poisson distribution, which is the case for count datasets such as term-document matrices or images. Most KL-NMF algorithms in the literature minimize a separable majorant of the loss to find their next iterate. We argue that this method has reached its limits and propose to use instead the second-order Taylor expansion of the loss, leading to a Newton-type method. We minimize this non-separable surrogate by proposing a generalization of the well-known HALS algorithm. This yields an efficient KL-NMF algorithm which provably converges and which competes favorably with state-of-the-art algorithms on a large variety of datasets.
We study exact Kullback--Leibler (KL) projection for low-rank factorizations whose two nonnegative factors have prescribed row marginals and a shared, learned column marginal. For arbitrary positive row marginals of equal total mass, the joint KL projection reduces exactly to a strictly convex gauge-fixed dual with only $r-1$ effective variables; its Hessian is a sum of categorical covariance terms and admits $O((n+m)r)$ matrix-free Hessian--vector products. The projection theorem is objective-independent. We then specialize this geometry to doubly stochastic (DS) graph learning through $W=U\operatorname{Diag}(g)^{-1}V^\top$, where row-simplex factors with a common column mass induce an exactly DS graph without materializing an $n\times n$ optimization variable. Combined with observed-edge sparse fitting, a stochastic anchor-reduced manifold regularizer, and Bregman backtracking, the resulting mirror-descent method preserves exact feasibility at every accepted step. Under a nonvanishing latent-mass condition, it satisfies sufficient decrease and an $O(1/N)$ mirror-stationarity bound, while strictly positive accumulation points are KKT stationary. Matched clustering experiments show competitive accuracy, feasibility residuals near numerical precision, and favorable anytime behavior without a dense learned graph.
Abstract.
This paper gives a necessary and sufficient condition for a nonnegative matrix that has an orthogonal nonnegative matrix factorization (ONMF) via characterization of the null space. We propose an optimization model to minimize the Frobenius norm of the product of a given nonnegative matrix and a variable matrix subject to the constraints defined by the necessary and sufficient condition. Moreover, we present an augmented Lagrangian algorithm for solving this minimization model and prove the global convergence to a stationary point. Two factor matrices for the ONMF of the given matrix can be easily obtained by the outputs of the algorithm. Preliminary numerical results using synthetic and real-world data with applications in clustering show that our approach outperforms some existing ONMF methods regarding accuracy and robustness.
Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as recognition that the authors have followed reproducibility principles valued by SIMAX and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/Qilun-Luo/ONMF . [Formula: see text]
High-dimensional multi-way (tensor) data pose significant challenges for covariance estimation due to the curse of dimensionality. We introduce a unified framework for scalable estimation of tensor covariances based on a Kronecker-structured sparse inverse Cholesky (KSIC) projection. Our approach is grounded in the geometry of information projection, defining the estimator as the moment-matching projection of a target distribution onto a manifold characterized by sparse, Kronecker-factored inverse Cholesky factors. By leveraging physical or data-driven nearest-neighbor sparsity, KSIC provides a geometry-aware representation that is both statistically interpretable and computationally efficient. Our framework integrates two estimation regimes: a nonparametric estimator that projects the empirical covariance directly onto the manifold, utilizing the KSIC structure to implicitly regularize rank-deficient data; and a parametric estimator that fits generative covariance models (e.g., Mat\'ern) by maximizing the likelihood of their KSIC projections, formulated as a nested double forward Kullback-Leibler minimization. Theoretically, we establish the conditions for the existence of the KSIC projection and finite-sample concentration rates for the nonparametric regime, proving that the KSIC estimator gainfully exploits cross-mode information and is robust to data scarcity. Numerical experiments demonstrate that the proposed KSIC estimators achieve state-of-the-art accuracy and scalability, particularly in settings with high dimensionality and limited sample sizes. We apply KSIC to spatiotemporal temperature anomalies and functional MRI data, demonstrating its broad applicability across diverse multi-way data domains.
We describe a simple rejection-sampling-based algorithm to perform length-squared sampling on an $n \times n$ positive-semidefinite (psd) matrix: that is, to sample a column with probability proportional to its squared $\ell_2$-norm. The algorithm runs in just $O(n)$ expected time, which is significantly sublinear in the input matrix size. The runtime is optimal, even when the input is assumed to be diagonal. Our result has several applications. Length-squared sampling is used by a number of sublinear time algorithms for matrix problems, like low-rank approximation and eigenvalue approximation. Often, it is assumed that the algorithm is given access to the matrix column norms, and thus can perform length-squared sampling efficiently. Our result shows that, at least for psd matrices, we can remove this assumption. We also discuss an application to an asymptotically optimal algorithm for estimating the Frobenius norm of a psd matrix to relative error. Finally, we show that our sampling algorithm yields a very simple sublinear time algorithm for the robust psd low-rank approximation problem introduced by Bakshi et al. (FOCS, 2020), which nearly matches the more complex method developed there.
Rajarshi Bhattacharjee, Ethan N. Epperly, Cameron Musco et al.· 0 citations
Laplacian-regularized minimization is fundamental in signal processing and machine learning, but is limited by the dense and ill-conditioned nature of the graph Laplacian pseudoinverse. While the Laplacian itself is sparse, its pseudoinverse is dense and often ill-conditioned, rendering direct computation impractical at scale. Moreover, pseudoinverse learning is more challenging than Laplacian learning. To address this challenge, this paper considers the setting where the graph Laplacian is given and proposes a Difference-of-Convex Regularizer (DCR) graph learning framework that approximates the spectral action of the Laplacian pseudoinverse without direct inversion via regularized Maximum Likelihood Estimation (MLE). By reformulating Laplacian-Regularized Nonnegative Least Squares (LR-NNLS) through a dual representation, DCR decouples pseudoinverse learning from instance-specific inference and enables efficient primal solution reconstruction via a differentiable dual-guided learning scheme. We establish theoretical guarantees on stability and the existence of a unique fixed point for DCR algorithm. Numerical experiments demonstrate improved performance over convex solvers and graph filtering baselines and robust performance across diverse graph topologies.
Nonnegative matrix factorization represents nonnegative signals as additive combinations of latent components, but its factorization rank, and hence the model order, must usually be specified beforehand. An underestimated order discards signal structure, whereas an overestimated order produces redundant components and unstable decompositions. We propose a column $\ell_{2,0}$-regularized formulation that estimates the model order from an initial upper bound by suppressing inactive columns in both factors. A warm-started regularization path progressively removes redundant components without changing the factor dimensions, and a marginal reconstruction-loss criterion selects an order along the path. To solve the resulting nonconvex and discontinuous problem, we develop an inertial proximal alternating linearized minimization method, a scale-balanced variant, and a proximal active-set method based on P-stationarity. The balancing operation equalizes the norms of paired factor columns while preserving their rank-one products. We characterize the critical points and local minimizers of the model, provide a sufficient-condition result for rank recovery, and prove whole-sequence convergence of the proposed algorithms under explicit step-size and inertial-parameter conditions using the Kurdyka--\L{}ojasiewicz framework. Dedicated experiments show that the warm-started $\lambda$-path is more efficient than increasing- and decreasing-order discrete $r$-paths, while scale balancing yields a more stable rank-selection path. Experiments on synthetic data and diverse signal benchmarks show that both iPALM and PASM provide reliable model-order estimates with favorable computational efficiency.