Fast randomized numerical rank estimation for numerically low-rank matrices

TL;DR

This work develops a fast, randomized approach to estimate the numerical rank of numerically low-rank matrices by sketching both column and row spaces to form $YAX$, enabling reliable tracking of the leading singular values $\sigma_i(A)$ through $\sigma_i(AX)$. It provides rigorous bounds showing leading singular values are preserved up to a constant factor under sketches, extends to general subspace embeddings, and analyzes tail-spectral effects and oversampling. The authors propose a practical numerical rank estimation scheme that uses a two-sketch pipeline to obtain $\hat\sigma_i$, integrates seamlessly with randomized low-rank approximation, and delivers speedups for fixed-precision problems by adjusting sketch sizes on the fly. Numerical experiments demonstrate robust gap detection, favorable comparisons to existing rank-estimation methods, and significant speed advantages in fixed-precision contexts, underscoring the method's value as a preprocessing step in large-scale low-rank computations.

Abstract

Matrices with low-rank structure are ubiquitous in scientific computing. Choosing an appropriate rank is a key step in many computational algorithms that exploit low-rank structure. However, estimating the rank has been done largely in an ad-hoc fashion in large-scale settings. In this work we develop a randomized algorithm for estimating the numerical rank of a (numerically low-rank) matrix. The algorithm is based on sketching the matrix with random matrices from both left and right; the key fact is that with high probability, the sketches preserve the orders of magnitude of the leading singular values. We prove a result on the accuracy of the sketched singular values and show that gaps in the spectrum are detected. For an $m\times n$ $(m\geq n)$ matrix of numerical rank $r$, the algorithm runs with complexity $O(mn\log n+r^3)$, or less for structured matrices. The steps in the algorithm are required as a part of many low-rank algorithms, so the additional work required to estimate the rank can be even smaller in practice. Numerical experiments illustrate the speed and robustness of our rank estimator.

Figures (9)

• Figure 1: How the tail of the singular value spectrum influences the singular value estimates. Each of the graphs shows the first 30 singular values $\sigma_i(A)$ of a square matrix of dimension $n=1000$. The first 20 singular values are identical for each of the matrices. The matrices in the graphs have constant, slow polynomially decreasing exponentially decreasing tails (from left to right). We sketch with a Gaussian embedding matrix of size $r= 19$. Although $\sigma_{r+1}=\sigma_{20}$ is the same for each of the matrices, we see the tail affects (only) the last few estimates. The results shown here are the average of 1000 trials, yet the behaviour is very typical.
• Figure 2: One could argue that in none of the cases discussed in Figure \ref{['fig:effecttails1']} is the numerical rank very well-defined, because there is no clear gap between any of the singular values. We repeat the experiment of Figure \ref{['fig:effecttails1']} but with slightly different singular values. The top row shows the estimates of the first 19 singular values when there is a gap between the 20th and 21st singular value of $A$. We see that the effects of different tails becomes negligible. In the bottom row we repeat the experiment with $r=25$ instead of $r=19$. This illustrates that we can detect the gap in each of the cases.
• Figure 3: Numerical rank estimation of matrices with polynomially decaying singular values and exponentially decaying singular values, where only the sketch $AX\in\mathop{\mathrm{\mathbb{R}}}\nolimits^{m\times r}$ is used. The matrices are real and square diagonal of dimension $n = 10^5$. The results are shown when using a subsampled randomized DCT embedding (SRTT), a Gaussian embedding matrix (Gaussian) and a hashed randomized DCT embedding (HRTT). The black lines indicate the exact numerical rank, $\text{rank}_{\epsilon}(A)$.
• Figure 4: Relative error in the singular value estimates of a real matrix $B$ with fast polynomially decaying singular values of size $10^5\times 2000$ based on the singular values of $Y B$, where $Y\in\mathop{\mathrm{\mathbb{R}}}\nolimits^{4000\times 10^5}$ is an embedding matrix (either a subsampled randomized DCT (SRTT), a Gaussian embedding matrix or a hashed randomided DCT (HRTT)).
• Figure 5: On the left, the singular values of $A_G$ and $YA_{G}X$, where $X$ is a Gaussian embedding and $Y$ is an SRTT matrix with a discrete cosine transform (SRCT). The plots from top to bottom correspond to $r_1 = 110, 210,310,410$ respectively. The gaps that exist in the spectrum of $A$ are also visible in the spectrum of $YA_{G}X$, even for values of $r_1$ close to the location of the gaps. On the right, the density of state plots resulting from the Chebyshev-based ubaru2016fast and Lanczos-based algorithms Ubaru2017. The algorithms are unable to identify the gaps for very small singular values. Even the exact regularized DOS can only identify the bulk of eigenvalues is of order $10^{-2}$ or smaller. See the references for further details on interpretation of these graphs.

Theorems & Definitions (14)

• Definition 1
• Lemma 2
• proof
• Theorem 3: Marčenko and Pastur Marcenko1967DistributionMatrices, Davidson and Szarek davidson2001local, Aubrun and Szarek aubrun2017alice
• Theorem 4
• proof
• Theorem 5
• proof
• Corollary 6
• Lemma 7