Across test matrices with diverse singular-value decays, the proposed methods attain ranks close to the truncated-SVD optimum while meeting the prescribed tolerance with high probability.
Abstract
We study fixed-tolerance low-rank approximation in the matrix-free setting, where a matrix or linear operator $\mathbf{A}$ is accessible only through matrix-vector products and its rank must be determined adaptively to meet a prescribed error tolerance. We introduce a family of adaptive, matrix-free randomized QB algorithms. A randomized error indicator estimates the residual norm -- in either the Frobenius or the spectral norm -- directly from a random sketch, remaining accurate down to machine precision. A matrix-free rank-pruning step decouples the computational block size from the final rank, so that large, BLAS-3-friendly blocks can be used without over-estimating the rank, and an adjoint-free variant returns the orthonormal basis using only the forward operator. Across test matrices with diverse singular-value decays, the proposed methods attain ranks close to the truncated-SVD optimum while meeting the prescribed tolerance with high probability.
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
We investigate when a sparse nonnegative matrix can be recovered from a real-valued matrix of much lower rank by zeroing out its negative elements. The potential for such decompositions suggests a mathematical connection between sparsity and rank; we analyze a number of sparse matrices with this latent low-rank structure and use them to illustrate the geometric origins of this connection. Previous algorithms have discovered these decompositions via an alternating minimization over the factors of a low-rank matrix, but to do so, they have also needed to compute and store another matrix, neither sparse nor low-rank, that is the size of their product. We develop a stochastic, alternating least-squares algorithm that operates on smaller blocks of this dense matrix and scales as a result to much larger problems. We also show how to further accelerate this algorithm with sparse optimizations and customized CUDA kernels. As one example, we use the algorithm to analyze the sparse matrix of synaptic weights for the recently published $\textit{Drosphilia}$ connectome. The nonzero elements of this matrix, with 139,255 rows and columns, record the number of synapses between cells in the nervous system of a female fruit fly. Despite a slowly decaying spectrum of singular values, this matrix exhibits a latent low-rank structure that is predictive of cell categories across multiple levels of specificity.
Lawrence K. Saul, N. Huang, Dennis Bollweg et al.· 0 citations
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.
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.
Quadratic optimization becomes hard as soon as either the matrix in the quadratic form has an unfavorable curvature or the feasible set is discrete, combinatorial, or otherwise nonconvex. A complementary phenomenon is also well known in the signal-processing and optimization communities: when the matrix in the quadratic form has small rank, some hard-looking quadratic programs admit exact polynomial-time algorithms for fixed rank. We study the common positive-semidefinite geometry behind this phenomenon. If $Q=BB^\top$ is positive semidefinite, the objective depends on $x$ only through the rank-space shadow $y=B^\top x$. Every optimal shadow $y^*$ uniquely maximizes the linear functional defined by its own direction and satisfies a quantitative quadratic margin. Thus nonlinear optimality collapses to a low-dimensional, self-generated linear exposure direction. We call this the rank-collapse principle. The principle alone does not imply a finite candidate set: efficient exact optimization additionally depends on the projected or active geometry of the feasible family. We organize this distinction through projected-shadow scattering and active-structure collapse, relate it explicitly to established zonotope, convex-combinatorial, edge-skeleton, projected-normal-fan, and fixed-rank sparse-PCA methods, and derive tie-safe consequences for binary and finite-phase vectors, cardinality constraints, matroid bases, and sparse PCA. We also give directional-stability and approximately low-rank certificates, together with reproducible experiments. The paper's contribution is a unified, careful framework and a set of quantitative consequences, rather than a claim to originate the known fixed-rank tractability results that motivate it.
Cubic regularized Newton methods have the optimal $\mathcal{O}(\epsilon^{-3/2})$ global rate, but a dense subproblem solve limits the feasible block size. Scalable Cubic Newton variants replace the true block curvature with a diagonal, low-rank, Kronecker-factored, or sketched surrogate and, most often, give up the exact cubic step. We introduce a blockwise optimizer that minimizes an independent cubic model per parameter tensor over the true block Hessian, under a per-block adaptive cubic constant and a monotone guard on the full loss. Arbitrarily large tensors are handled matrix-free in a Lanczos-built Krylov subspace, where we prove that the step minimizes the cubic model. The theory also supplies the $\mathcal{O}(\epsilon^{-3/2})$ iteration complexity bound, a second-order guarantee, and monotone per-block descent. Four variants of this outer scheme are evaluated against the original adaptive regularization with cubics (ARC) optimizer, some other recent cubic Newton variants, Adam, SOAP, and L-BFGS. On a 91.4M-parameter implicit neural representation (INR), the variants introduced in this work are the only evaluated here cubic Newton methods whose steps stay exact on every block. Run to full convergence on FINER 2D image fitting, one of the ARC variants introduced here, ARC-$\varphi_1$, reaches 133.5 dB peak signal-to-noise ratio, while tuned Adam plateaus at 78.2 dB after about 70 minutes. In that time ARC-$\varphi_1$ reaches 95.6 dB.