Uniform Sampling from High-dimensional Spectral Norm Balls

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.

Non-Euclidean optimisation methods with matrix-valued updates, such as Muon and Scion, have recently shown strong empirical performance for training Transformer models, yet their theoretical advantages over Euclidean methods remain poorly understood. We address this gap in the heavy-tailed non-convex regime, where stochastic gradients have bounded $p$-th central moments, $p \in (1,2]$. We show that certain non-Euclidean methods achieve optimal sample complexity under stronger stationarity measures, while Euclidean methods incur additional dimension-dependent costs. As a consequence, for $m \times n$ matrices, Muon finds an $\varepsilon$-stationary point in nuclear norm within $\mathcal{O}\left(\min\{m, n\} \frac{Δ_1 L}{\varepsilon^2} \left(\frac σ\varepsilon \right)^{\frac p {p-1}}\right)$ samples, absorbing heavy-tailed noise without extra dimension dependence, unlike Euclidean methods. We further prove this sample complexity, including its dimension dependence, is optimal for all first-order methods under nuclear-norm stationarity. Experiments on large language models support our theory. Surprisingly, our results suggest that other Schatten geometries beyond the spectral geometry of Muon can perform competitively in certain settings.

Over the past half-dozen years, stochastic rounding (SR) has regained significant attention as a quantization scheme for low-precision floating-point arithmetic, with applications spanning numerical analysis and modern machine learning systems. Recent work has shown that SR acts as an implicit regularizer by increasing the smallest singular value of extremely tall-and-thin (or, symmetrically, short-and-fat) matrices. In this work, we substantially sharpen and extend this understanding in two directions. First, we show that the regularization effect of SR is not restricted to extreme aspect ratio regimes: it persists for matrices with constant aspect ratio. Second, we demonstrate that SR does not merely regularize the smallest singular value, but instead lifts entire clusters of singular values at the tail of the spectrum. Together, these results provide a more general characterization of stochastic rounding as a spectral regularizer, revealing that its effects extend beyond extremal aspect ratios and act on a broader portion of the singular value spectrum.