Fast randomized numerical rank estimation for numerically low-rank matrices

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.


Introduction
Low-rank matrices are ubiquitous in scientific computing and data science.Sometimes a matrix of interest can be shown to be of numerically low rank [6,45], for example by showing that the singular values decay rapidly.More often, matrices that arise in applications may have a hidden low-rank structure, such as low-rank off-diagonal blocks [30,32].As is well known, low-rank approximation is also the basis for principal component analysis.
Numerous studies and computational algorithms exploit the (approximately) low-rank structure of these types of matrices to devise efficient algorithms.Such algorithms usually require finding a low-rank approximation of a given matrix.For large-scale matrices, traditional deterministic algorithms for low-rank approximation based on computing a truncated singular value decomposition (SVD) might be infeasible.This can be caused by the sheer scale of the matrix and the complexity of classical algorithms, as well as the fact that access to the matrix may be restricted because it is stored in RAM or cannot be stored at all (streaming model) [32].Randomized algorithms have become a powerful and reliable tool for efficiently computing a near-best low-rank approximation for such matrices.Landmark references on randomized low-rank factorizations are [48] and [19], which extensively analyze the randomized SVD.
A key step in low-rank algorithms, including randomized SVD, is the determination of the numerical rank (based on the spectral norm).For instance, a variety of low-rank factorization techniques require the target rank of the factorization as input information.Once a low-rank approximation Âr of the specified target rank has been computed, a posteriori estimation of the error ∥A − Âr ∥ 2 via a small number of matrix-vector multiplications [32] is a reliable means of checking if the input rank r was sufficient.If r was too low, one would need to sample the matrix with more vectors.This clearly requires more computational work, and can be a major difficulty in the streaming model, wherein revisiting the matrix is impossible [41].Conversely, if the input rank r was too high, the computational cost of computing Âr would be higher than it could be.Having a fast and reliable rank estimator is therefore highly desirable.
In this work we propose and analyze a fast randomized algorithm for estimating the numerical rank of an approximately low-rank matrix A ∈ R m×n or A ∈ C m×n .The algorithm is based on sketching both the column as well as the row space of A, forming Y AX, where Y and X are random (subspace embedding) matrices.The key idea is that with high probability, the singular values of AX are good estimates of the (leading) singular values σ i (A) of A. Therefore the decay of σ i (A) can be reliably estimated by the decay of σ i (AX).To estimate σ i (AX) we once again sketch AX to obtain the much smaller matrix Y AX; it is only σ i (Y AX) (or their estimates) that we actually need to compute.It is noteworthy that in the algorithm it is only necessary to view to matrix A once.Furthermore, since we obtain estimates σ i (Y AX) for the leading singular values σ i (A), we will be able to detect a gap in the spectrum of A and hence have the possibility of setting a tolerance for the numerical rank in a data-driven manner.
To our knowledge, the fact that a sketch of the form Y AX preserves 'some coarse spectral structure' was first noted in [1].We use techniques similar to their proofs for our theoretical results.Related results for the bulk of the eigenvalues of (AX) * (AX) were first obtained in [17].The main theoretical contributions of this work, deterministic and probabilistic bounds on the accuracy of σ i (AX) as estimates for σ i (A), build on these results.In particular, we provide a tighter lower bound on σ i (AX)/σ i (A) relevant for small singular values of AX.
The complexity of our algorithm is O(mn log n+r 3 ) for dense m×n matrices, and can be lower if A has structure that can be exploited for computing the sketches AX and Y AX.This is clearly a subcubic complexity (assuming r ≪ m, n), and it runs significantly faster than computing the full SVD.Moreover, in many cases, computing AX (which is usually the dominant part of our rank estimation algorithm) is a required part of the main algorithm (e.g.randomized SVD); and in some algorithms [11,36] this is true even of Y AX.In such cases, the additional work needed for estimating the rank is therefore significantly smaller, such as O(r 3 ) or sometimes even O(r).Our algorithm is competitive when the target rank r is much smaller than the dimensions of the matrix.
In particular, for the widely used randomized SVD we argue that it is computationally efficient to compute Y AX and its singular values, after having computed a sketch AX.The reason is that σ i (Y AX) will provide information on how appropriate the size of the sketch AX is.Subsequently, more columns can be added or columns can be removed before computing a QR decomposition QR = AX and computing the product B = Q * A. This will save computational cost in either case: an unnecessary large sketch, or a sketch of insufficiently large dimension.We also introduce an algorithm that combines the rank estimation algorithm with randomized SVD, using the information provided by σ i (Y AX).
Notation.Unless specified otherwise, A is an m × n matrix where m ≥ n.The ith largest singular value of A is denoted by σ i (A), and we furthermore use σ max (A) and σ min (A) for the largest and smallest singular values respectively.We use ∥•∥ to denote the spectral norm ∥A∥ = σ 1 (A), and ∥A∥ F is the Frobenius norm.F denotes the field, in our case F = R or F = C.The numerical rank estimate will be denoted by r.
Throughout the paper we use X and Y for random oblivious subspace embedding matrices, such that with high probability, for any fixed Q with orthonormal columns we have σ max (Q * X), σ max (Y Q) ≤ 1 + ε and σ min (Q * X), σ min (Y Q) ≥ 1 − ε for some ε ∈ (0, 1).We use r, r 1 and r 2 to refer to the size of embedding matrices, which must be at least the number of columns in Q.The matrix X is required to have more columns than Q, and the same goes for the rows of Y .For brevity we simply call such X and Y embedding matrices.
The analysis will be specified to Gaussian or subsampled randomized trigonometric transform (SRTT) [32] embedding matrices at times.We will use G ∈ R n×r , where n > r, to denote a matrix where each entry is iid N (0, 1) distributed and call this matrix a real (standard) Gaussian matrix.In case of F = C, G ∈ C n×r denotes a complex Gaussian matrix with iid CN (0, 1) entries.That is, Re(G ij ), Im(G ij ) ∼ N (0, 1/2) and independent.A Gaussian matrix can be scaled to an embedding matrix by defining X = G/ √ r (in both the real and complex case); we call the scaled Gaussian matrices Gaussian embedding matrices.SRTT matrices are defined in Section 3.3 and as the scaling is included in the definition, they are naturally embedding matrices.We will use Θ to denote an SRTT matrix.

Numerical rank and goal of a rank estimator
So far we have been using the term "the numerical rank" informally.Let us now define the notion.This is a standard definition, see for example [6].
We adopt the relative ϵ-rank definition as it is a natural choice in the context of low-rank matrix approximation.One could easily adapt the algorithm presented to estimate the absolute ϵ-rank.Let us discuss what the goal of a rank estimator should be.One natural answer of course is that the estimator should output rank ϵ (A) given A and the (user-defined) relative threshold ϵ.However, we argue that the situation is more benign: consider, for example, a matrix with ∥A∥ = 1 and k singular values > 10ϵ, five singular values in (ϵ, 1.1ϵ), five more in (0.9ϵ, ϵ), and the remaining n − k − 10 are < 0.01ϵ.What should the estimator output?Is it crucial that the "correct" value rank ϵ (A) = k + 5 is identified?
In this paper we take the view that the goal of the rank estimator is to find a good ϵ-rank, not necessary the ϵ-rank.In the example above, any number between k and k + 10 would be an acceptable output.In most applications that we are aware of (including low-rank approximation, regularized linear systems, etc), there is little to no harm in choosing a rank r ̸ = rank ϵ : a slight overestimate r > rank ϵ usually results in slightly more computational work, whereas a slight underestimate r < rank ϵ is fine if σ r+1 (A) = O(ϵ∥A∥), that is, a rank-r matrix can approximate A to relative O(ϵ)-accuracy.
On the other hand, it is clearly not acceptable if the rank is unduly underestimated in that σ r+1 (A) ≫ ϵ∥A∥.It is also unacceptable if σ r (A) ≪ ϵ∥A∥.The goal in this paper is to devise an algorithm that efficiently and reliably finds an r such that • σ r+1 (A) = O(ϵ∥A∥) (say, σ r+1 (A) < 10ϵ∥A∥): r is not a severe underestimate, and • σ r (A) = Ω(ϵ∥A∥) (say, σ r (A) > 0.1ϵ∥A∥): r is not a severe overestimate.
Any such r is sensible in that there exists a rank-r approximation of A with O(ϵ) relative accuracy, and r is not much larger than the best possible for the accuracy required.In addition, the nature of randomized algorithms means that the user must be willing to accept results that do not hold precisely and that small inaccuracies in the results are acceptable.As a result, we assume one is not looking for the exact numerical rank, but an order-of-magnitude estimate.Our rank estimate will satisfy these conditions with high probability.In particular, in situations where the numerical rank is clearly defined, i.e., if a clear gap is present in the spectrum σ r (A) ≫ ϵ∥A∥ ≫ σ r+1 (A), the algorithm will reliably find the exact rank r = rank ϵ (A).
Additionally, in situations where it is unknown whether the matrix A is lowrank approximable, our rank estimator can tell us (roughly) how well A can be approximated with a low-rank matrix.Similarly, if it is unknown what a suitable tolerance would be and/or the user would like to detect a gap in the spectrum, the algorithm can be used to plot an estimated spectrum and detect gaps.
This paper is organized as follows.Section 2 discusses related studies in the literature.In Section 3 we show that σ i (AX) gives useful information about σ i (A) for leading values of i.Then in Section 4 we show that σ i (AX) can be estimated via σ i (Y AX).Section 5 summarizes the overall rank estimation algorithm.Numerical experiments are presented in Section 6.

Existing methods for numerical rank estimation
At the core of numerical rank estimation lies singular value estimation; specifically, the rank can be found by counting the number of singular values greater than the tolerance ϵ.The most direct way to do this is to explicitly compute the singular values of a matrix.However, as is well known, this requires O(mn 2 ) operations for an m × n matrix [18], and for large matrices this is computationally inadmissible.Another point of view, specifically for Hermitian matrices, is that a numerical rank estimate can follow from estimating the density of states (DOS).The DOS can be interpreted as a probability density function describing the position of the eigenvalues, and can be estimated with algorithms analyzed in [27].
Regarding fast algorithms that can run with subcubic complexity, the literature on rank estimation appears rather scant.Exceptions include the work of Ubaru and Saad [42], Ubaru, Saad and Seghouane [43], Zhang, Wainwright and Jordan [50] and Di Napoli, Polizzi and Saad [14].These works are all concerned with counting the number of eigenvalues in a certain interval for a Hermitian matrix.Alternatively, [10] discusses an algorithm for computing the exact rank of a matrix.
The first two references [42,43] employ the idea of density of states (DOS).In [42], the authors use the DOS to locate a gap in the spectrum, derive an appropriate tolerance ϵ and subsequently use polynomial approximation and stochastic trace estimation to count the number of eigenvalues greater than the tolerance.Computationally, the algorithm only requires matrix-vector products with A but consequently also requires many views of A. In the second paper [43], the authors also first estimate the DOS to locate a gap and then estimate the integral of the DOS.We compare against these algorithms in the numerical experiments, although they need to be adjusted slightly to account for matrices that are not symmetric positive semi-definite.An advantage of both methods is that the entire spectrum is estimated, however they require many more views of the matrix than most randomized methods.The results are discussed in Section 6.3.The works [50] and [14] are similar in spirit.The first focuses on minimizing communication cost involved with rank estimation in a distributed setting, the second approximates an eigenvalue count using polynomial and rational approximation filtering.
As for more theoretical contributions, much of this work builds upon Andoni and Nguy ễn [1], which is to the best of the authors' knowledge the first work that shows rank can be estimated via sketches.The work is focused on estimating eigenvalues of Hermitian matrices via the sketch X T AX, but also considers the singular values of non-symmetric matrices using Y AX.While Andoni and Nugyen work with the Gram matrix of the sketch, obtaining results in terms of the difference of the squared singular values (or via the Jordan-Wieldant matrix, which is larger), we work with Y AX directly.Moreover, their work is of a theoretical nature, and do not present a concrete algorithm for estimating the rank.They also state results that only come with relatively low probability guarantees, although it is possible to take a larger sample to improve the probability.
Here we introduce precise algorithms (also for low-rank approximation), show numerical experiments, and phrase our results in terms of popular embedding matrices allowing for results that hold with exponentially small failure probability.

Related work in statistics
Rank estimation, or singular value estimation, is a well-known problem in the context of covariance matrix estimation or PCA in statistics, or the detection of signals in signal processing.Most rank estimation algorithms in this context require computing an SVD or eigendecomposition; exceptions include [44].Since we consider a context where computing these decompositions exactly is not computationally feasible, we do not take this into account.
The statistics literature is connected to this work in yet another way: relating the singular values of the sketch AX to the singular values of A is equivalent to the problem of covariance matrix estimation.This can be seen as follows: let A = U ΣV * be the SVD of A and note the singular values of AX are the same as the singular values of ΣX 1 , where X 1 = V * X is a Gaussian embedding matrix.The columns of ΣX 1 can be viewed as r scaled observations from an n-dimensional N (0, Σ 2 ) distribution.As the number of observations tends to infinity, the singular values of the 'data matrix' ΣX 1 tend to the singular values of Σ. Estimation of (approximately) low-rank covariance matrix is studied in, among others, [46,47,23] and has a clear relation to PCA.
One specific related example in the statistics literature that is well studied is the spiked covariance matrix model, introduced by Johnstone [21] and analyzed by, among others, Bai and Silverstein [4], Nadler [35] and Rao et al. [37].We find that the theory in these works cannot be applied in our context, since the tail of the singular values in this model is very heavy.In the general context of covariance matrix estimation, various authors have suggested different manipulations to the sample covariance matrix or its eigenvalues to improve the estimator, known as shrinkage.See, for instance, the work by Ledoit and Wolf on linear shrinkage [24], nonlinear computational shrinkage [25] and nonlinear analytical shrinkage [26], or by El Karoui based on random matrix theory [16].
Our numerical experiments suggest these methods are either inefficient for largescale matrices or unsuitable for numerically low-rank matrices, and we do not discuss them further.

Related algorithms in (randomized) numerical linear algebra
A number of recent papers focus on randomized algorithms to perform numerical linear algebra tasks on large matrices.Specifically, randomized algorithms for either full or low-rank matrix factorization are of interest in the discussion of numerical rank estimation.This relevance is two-fold, as a numerical rank estimate can be a natural consequence of a factorization.In other cases, an a priori rank estimate can aid in performing a (low-rank) factorization.
The former is the case for rank-revealing full factorizations, such as rankrevealing QR decompositions [15] or UTV factorizations [31].Full factorizations are (relatively speaking) most useful when the ϵ-rank is not much smaller than the matrix dimensions, and still require cubic O(mn 2 ) complexity.As we mainly consider the case where rank ϵ (A) ≪ min(m, n), we focus on methods that have subcubic complexity.
A number of randomized algorithms for a low-rank factorization have been proposed and they can often be modified to yield a (rank r) QB factorization: A ≈ QB, where Q ∈ F m×r is a matrix with orthonormal columns and B ∈ F r×n is a short and fat matrix.Singular value estimates can then be obtained by computing the exact singular values of B, which in turn leads to a numerical rank estimate.Within randomized algorithms for QB factorizations, we can distinguish between (adaptive) algorithms that focus on solving the fixed-precision problem and non-adaptive algorithms that focus on the fixed-rank problem.
The fixed-precision problem can be summarised as: find QB such that ∥A − QB∥ ≤ ϵ∥A∥.Adaptive algorithms for this problem include the adaptive randomized range finder [19], the incremental rangefinder [32], RandQB_b [29] and RandQB_EI [49].The fixed-precision problem is discussed in more detail in Section 6.4.
Conversely, the fixed-rank problem is: find QB of exact rank r such that ∥A − QB∥ is as small as possible.One of the most well-known algorithms in randomized fixed-rank factorization is randomized rangefinder [19], which is a core part of the randomized SVD algorithm proposed in [19].
We discuss these algorithms further, comparing with our algorithm, in Sections 5-6, to demonstrate that our algorithm is much more efficient for rank estimation, and can often be used as a convenient preprocessing step for a lowrank approximation algorithm to determine the input rank.
Conversely, a numerical rank estimate can aid the factorization process.This is discussed in detail in Section 6.4 and forms the basis for a non-adaptive lowrank factorization algorithm we propose for the fixed-precision problem.

Sketching roughly preserves singular values: σ i (AX)/σ i (A) = O(1)
for leading i In this section we show that if a matrix A has low numerical rank, then sketching preserves its leading singular values sufficiently well to provide an estimate for its numerical rank.Again drawing the connection to sample covariance matrix estimation, it is intuitive that it should be difficult, if not impossible, to obtain information about an n-dimensional distribution using only r samples, when r ≪ n.Clearly, there are n − r singular values (or signals) we cannot detect, but it is also not clear why the r singular values we estimate would be any good.One way to think about this is to see that N (0, Σ 2 ) lies close to a low-dimensional subspace, since A is approximately low-rank.This makes it more intuitive why a small sample size could suffice to retrieve a reasonably accurate approximation of Σ (see also [46,Sec. 9.2.3]).
We make the connection between the low-rank structure of A and our ability to estimate its singular values using AX explicit with the deterministic and probabilistic error bounds on the ratio σ i (AX)/σ i (A).The analysis will focus on Gaussian matrices.As explained before, G ∈ R n×r denotes a Gaussian matrix with iid N (0, 1) (F = R) or CN (0, 1) (F = C) entries, and a Gaussian embedding matrix will be of the form X = G/ √ r.We start from the unscaled sketch AG.Let A = U ΣV * be the SVD of A and decompose this further to where U 1 contains the leading r singular vectors of A. Now decompose AG ∈ F m×r as where Σ 1 is r × r and G 1 ∈ F r×r and G 2 ∈ F (n−r)×r are independent Gaussian matrices.Define We start with the following result, which examines the relation between σ i (AG) and σ i (A).

Lemma 2.
Let A ∈ F m×n and G ∈ F n×r .Decompose AG as in (1).Then, for i = 1, . . ., r, where Ĝ{i} is an i × r matrix consisting of the first i rows of G 1 , and G{r−i+1} is the matrix of the last r − i + 1 rows of G 1 .Furthermore, if G is a Gaussian matrix, then for each i ∈ {1, . . ., r}, the random matrices Ĝ{i} , G{r−i+1} , and G 2 are independent Gaussian in F.
Proof.The proof consists of establishing the identities for i = 1, . . ., r.
It is immediate from the definitions of B 1 and As a result, we have by Weyl's theorem, for i = 1, . . ., r, The upper bound follows from observing The lower bound is immediate.We use the interlacing property of singular values to show (4).For a square matrix B ∈ F k×k , let B [−p] ∈ F k×(k−p) denote a submatrix of B obtained by deleting any p columns.Then [33] First, notice the upper bound ( 4) is trivial for i = 1.For i = 2, . . ., r, partition where Σ 11 ∈ R (i−1)×(i−1) and Σ 12 ∈ R (r−i+1)×(r−i+1) are diagonal matrices which contain the singular values of Σ 1 .We then have where we applied (6) with j = 1 and p = i − 1 for the first inequality.Similarly, for the lower bound and i = 1, . . ., r − 1, consider the partition where now Σ 11 is i × i and Σ 12 is (r − i) × (r − i).We have The lower bound is immediate for i = r.
If G is Gaussian, its orthogonal invariance implies G is also Gaussian.
Since Ĝ{i} , G{r−i+1} , and G 2 are submatrices of this matrix with no overlap for a fixed i, it immediately follows that they are Gaussian and independent.
It is worth noting that (2) holds without the assumption that G is Gaussian.Lemma 2 motivates us to look at the singular values of independent Gaussian matrices Ĝ{i} , G{r−i+1} and G 2 .To this end we use two classic results from random matrix theory to derive probabilistic bounds on the expected error and on the failure probability.
Furthermore, for every t ≥ 0 one has Let G ∈ C m×n be a complex Gaussian matrix with m ≥ n.Then, for every t ≥ 0 one has Furthermore, for every t > 4 √ 2 log n/( m/n − 1), A useful way to interpret these results in simple terms is that rectangular random (Gaussian) matrices with aspect ratio m/n > 1 are well-conditioned, with singular values supported in [ We note that (7), which follows from the Marčenko and Pastur rule [28], holds more generally for random matrices with iid entries with mean 0 and variance 1.The second pair of bounds (8) do not generally hold for other random matrices, see [39] for details.As far as the authors are aware, there does not exist a lower bound on the expectation of the smallest singular value of a complex normal matrix.An upper bound for the largest singular value can be found in [2,Prop. 6.33].We use these results to bound the singular values of Ĝ{i} , G{r−i+1} , and G 2 , which leads to the following theorem.We present the result only in the real case, but it can easily be extend to the complex case by combining Lemma 2 with ( 9) and (10).Theorem 4. Let X ∈ R n×r be a Gaussian embedding matrix, i.e., X = G/ √ r, and let A ∈ R m×n , where m ≥ n.Then, for i = 1, . . ., r Additionally, for each i = 1, . . ., r we have, with failure probability at most 3e −t 2 /2 , Proof.Note first from Marčenko-Pastur (7) that we have We can combine this with the result of Lemma 2 to find and The first result of Theorem 4 then follows by applying the scaling 1/ √ r.Similarly, by (8) we have with failure probability at most 3e −t 2 /2 that the following statements hold simultaneously: Since the deterministic results in Lemma 2 imply the result follows.

Bounds for general subspace embeddings
The results above focus on Gaussian embedding matrices, and show that the ratios of the singular values σ i (AX)/σ i (A) are reasonably close to 1, as specified by (11).We now show that much of this carries over to a general subspace embedding X.Note that an r-dimensional subspace embedding will not generally embed an r-dimensional space, but rather a smaller, say r-dimensional, space.Theorem 5. Let Ṽ1 ∈ F n×r be the leading r(≤ r) right singular vectors of A ∈ F m×n , and suppose X ∈ F n×r is a subspace embedding for the subspace Ṽ1 Proof.We start with a variant of (2), which follows from the same arguments: Here X{i} , X{r−i+1} are submatrices of Ṽ T 1 X (as in Lemma 2), and X 2 = Ṽ T 2 X, where Ṽ2 ∈ F n×(n−r) is the trailing right singular vectors of A. By the assump- which gives the lower bound in (13), and also to obtain the upper bound.We note that the ∥X∥ 2 term in ( 13) is bounded by n/r, often deterministically, in most commonly-used embeddings (see, e.g. ( 15)).This is roughly equal to n−r r appearing in (4).Thus Theorems 4 and 5 offer roughly the same level of guarantee in terms of the quality of σ i (AX) as an estimate for σ i (A).

The effect of tail singular values and oversampling
The proof of Lemma 2 allows us to see how the tail of the singular values, Σ 2 , affects the lower bound on the ratio σ i (AG)/σ i (A).This allows us to make a case for oversampling based on singular value estimation, whereas usually oversampling is inspired through a discussion on singular vectors.Again using the Courant-Fischer theorem, we have x .
Now let Ŝ ⊂ F r and x ∈ Ŝ be such that max We can use this to show As G 2 x is a Gaussian vector, we know E∥Σ 2 G 2 x∥ = ∥Σ 2 ∥ F .Additionally, it is worthwhile to note σ min ( Ĝ{i} )/ √ r ≈ 1 − i/r in expectation.The lower bound (14) thus suggests that a heavier tail Σ 2 , in terms of ∥Σ 2 ∥ F , would result in larger values for σ i (AX) than a light tail.Consider in particular the last few singular values of AX, where the term σ min ( Ĝ{i} ) is small.If the tail is light too, i.e. ∥Σ 2 G 2 x∥ is small, the lower bound will not be strong.We see in practice that this effect is indeed noticeable for the smallest singular values of AX, making them less reliable estimators.The bulk of the singular values of AX are, however, not affected by the size of the tail given the low-rank structure of the matrix.We display this effect in Figures 1 and 2 with smallscale experiments.
The small-scale experiments agree with bigger numerical results, as well as with theoretical results, in the sense that the extreme (smallest) singular value estimates are not trustworthy.For this reason we recommend oversampling by 10% and disregarding the singular values corresponding to the oversampling.

Choice of sketch matrix
In our analysis we have mainly focused on Gaussian embedding matrices, as they are the most well-studied class of random matrices and allow us to derive sharp constants.A Gaussian embedding matrix is a random matrix with iid  1 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 1 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.
entries N (0, 1/r).They are often the simplest sketch to implement.While it requires O(mnr) operations to compute AX, when r = O(1) the complexity is optimal as mn is a lower bound for dense A, as clearly all entries of A need to be read.In fact, Gaussian sketches can be among the fastest to execute when r = O(1).Furthermore, Gaussian is the most efficient type of embedding in the sense that an embedding for an r-dimensional subspace can be achieved with high probability using O(r) samples.
In particular, SRTTs such as the subsampled randomized Hadamard, Fourier, discrete cosine or discrete Hartley transforms are random embeddings that allow fast application to a matrix.The use of these sketches can be justified by Theorem 5.While they come with weaker theoretical guarantees than Gaussian embeddings [40] (in that the size of the sketch needs to be at least O(r log r) to ensure an embedding with high probability), empirical evidence strongly suggests that they usually perform just as well.The performance of SRTT matrices is discussed more in the next section.
Related to SRTTs are the recently introduced class of random embeddings, HRTTs [9].In these works it is shown an embedding for an r-dimensional subspace can be achieved with high probability using O(r) samples (as with Gaussians) rather than O(r log r), by replacing the subsampling matrix S in (15) with a hashing matrix (e.g.Countsketch [48]), which can be done without increasing the complexity.Computational results suggest HRTTs perform adequately in pathological cases where SRTTs fail, such as diagonal matrices (see discussions on coherence in [20,8]).
Finally, highly sparse embedding matrices have attracted much attention in the theoretical computer science literature [11,22].The state-of-the-art result [12] suggests an oversampling by a factor log r is needed to obtain an embedding with high probability.
When A ∈ F m×n has no structure to take advantage of (e.g.dense), such SRTT embeddings allow us to compute AX with O(mn log n) operations, which is optimal up to an O(log n) factor.When A has structure such as sparsity, AX could be computed much more efficiently, for example in O(nnz(A)r) operations for a Gaussian X and less for sparse sketches.
The choice of X therefore depends on the structure of A; but the complexity of the sketching AX can be bounded from above by O(mn log n) for dense matrices.We mainly use HRTT embedding matrices -specifically hashed randomized DCTs -in this paper as they are able to perform optimally, even for very coherent matrices.In our numerical experiments, we use diagonal matrices which is the most coherent type of matrix, and therefore most difficult to sketch.

A numerical experiment
We show how the theoretical result of this section, namely σ i (AX)/σ i (A) = O(1) for leading i, practically translates to a rank estimation technique.After sketching with an n × r embedding matrix X, we compute the singular values σ 1 (AX), . . ., σ r (AX) and estimate rank ϵ (A) for a given tolerance ϵ to be the first r such that σ r (AX) > ϵ∥A∥ ≥ σ r+1 (AX).If (an estimate of) ∥A∥ is unavailable, we may use σ 1 (AX).We applied this method to two synthetic example matrices, one diagonal matrix with polynomially decaying singular values (σ i (A) = i −3 ), and one diagonal matrix with exponentially decaying singular values (σ i (A) = 10 −0.01(i−1) ).The matrices are square of dimension n = 10 5 .
Figure 3 shows the results of the experiment.We plot the estimates for the numerical rank r that result from varying embedding dimensions r, and the ratio σ r+1 (A)/σ rankϵ(A)+1 .As both matrices do not have large gaps in the spectrum, the latter is an important quantity to judge the effectiveness of the rank estimator (see Section 1.1).We see that even for small values of r and when using an HRTT or Gaussian embedding, σ r+1 (A) is very close to σ rankϵ(A)+1 and ϵ and well within the acceptable range described in Section 1.Additionally, we see that SRTT matrices perform less than optimal for this pathological example of a diagonal matrix.

Randomized approximate orthogonalization: σ i (Y AX)/σ i (AX) = O(1)
We have seen that the numerical rank of AX is a good estimator for the numerical rank of A, provided that A is approximately low-rank.Although AX is of much smaller dimension than A, it may still be prohibitively expensive to compute its exact singular values.Furthermore, one could question whether it is worth calculating the singular values accurately spending O(mr 2 ) work, considering the accuracy that was lost in the sketching step.Here we describe a method, which we call randomized approximate orthogonalization, to cheaply compute estimates of the singular values of AX in O(mr log(m)+r 3 ) operations.The idea is similar to Section 3, in that sketches roughly preserve (the leading) singular values, but as we are working with a tall-skinny matrix AX, here all the singular values will be preserved up to a small multiplicative factor.
Randomized approximate orthogonalization is inspired by the randomized least-squares solver framework named sketch-to-precondition.The ideas were introduced in [38] and a fast implementation was introduced in Blendenpik [3].In the framework, a preconditioner for an overdetermined least-squares problem min x ∥Bx − b∥ where B ∈ R m×r1 (m ≫ r 1 is generated.The solvers work as follows: the first step is to sample the rowspace of B by sketching with a random embedding matrix Y ∈ R r2×m , where r 2 ≥ r 1 , to obtain the small matrix Y B. Secondly, one computes the QR factorization of Y B = QR.The upper triangular factor R is then used as a preconditioner for an iterative solver such as LSQR.Even though the QR factorization is based on a sketch of B, it is now a well-known fact in RNLA that BR −1 is 'close to orthonormal' in the sense that κ(BR −1 ) = O(1).Consequently, the singular values of R, or equivalently Y B, are close to the singular values of B (in the relative sense 1) for all i).One could also view this process as a 'randomized approximate orthogonalization' of B.
It follows that, in our context, we may estimate the singular values of AX using the singular values of Y AX, where Y is another random embeddding (independent of X) of size r 2 × m.It is then only necessary to compute the numerical rank (by computing the exact singular values) of the very small matrix Y AX ∈ F r2×r1 to estimate the numerical rank of the large matrix A. This is the final step of our rank estimation algorithm, which is discussed fully in the next section.
As for the choice of Y : unlike the first sketch of computing AX, in most cases the structure of A (such as sparsity, if present) is lost in B = AX.It is therefore usually recommended that we take Y to be an SRTT embedding, so that Y B = ΘB can be computed in O(mr log m) operations.We assume this choice below and throughout the remainder of the paper, switching between the notations Y AX and ΘAX when appropriate.That is, where S ∈ R r2×m is a sampling matrix, F a square orthogonal (or unitary) trigonometric transform (such as Fourier or DCT) of dimension m and D a diagonal matrix of independent random signs.Specifically the subsampled randomized Hadamard transform is analyzed extensively in [40,7].These results can be readily extended to a general SRTT matrix and used in our context, as we show in the appendix in Lemma 7.This is essentially a restatement of [40,Thm. 3.1] and [7,Lem. 4.1].It leads to the following result on the accuracy of σ i (Y AX) as estimates for σ i (AX).Corollary 6.Let AX ∈ R m×r1 , with m ≥ r 1 , and let Θ ∈ R r2×m be an SRTT matrix as in (15), where the trigonometric transform then, with failure probability at most 3δ The use of the number η = m max |F ij | 2 is a natural choice as it only depends on the choice of trigonometric transform, not on the size of the matrices involved.Choices for the trigonometric transform include Fourier, cosine, Hadamard, and Hartley transforms.The optimal η is attained for a Hadamard transform, for which η = 1.The same holds for a Fourier matrix, but this involves operations in C and when A is real, we may prefer real transforms.Since the Hadamard transform is only available when the dimension m is a power of 2 (for which one solution is to append zeros to the matrix), we make use of the discrete cosine transform with η = 2 in our experiments when A is real.Alternatively, we could have chosen the discrete Hartley transform with the same coherence η = 2. See [3] for more discussion.
We also note that a result analogous to Corollary 6 is given in [34,Thm. 4.4] when Θ is a Gaussian matrix (and hence in a closer spirit to Section 3).Results for a general subspace embedding matrix can be found in [48].

A numerical experiment
We explained how all singular values of a tall and skinny matrix AX ∈ F m×r1 can be estimated by the singular values of Y AX, where Y ∈ F r2×m is an embedding matrix and r 2 ≥ r 1 .In the next section we conclude how this, combined with the results of the previous section, leads to a rank estimation algorithm.In the following numerical experiment, we first illustrate the error one can expect in practice in the step from AX to Y AX.
We let B be a dense matrix with fast polynomially decaying singular values, σ i (B) = i −3 , of size 10 5 × 2000.We estimate σ i (B) with σ i (Y B), where Y is either a subsampled randomized DCT embedding or a Gaussian embedding of dimension 4000.The relative error |σ i (Y B) − σ i (B)|/σ i (B) is plotted in Figure 4. We see the relative errors are o(1) for each i = 1, . . ., 2000.We use a dense matrix with incoherent left singular vectors to resemble the matrix AX, which will never be very coherent due to the prior application of embedding X.As a result, the SRTT embedding is sufficient. is an embedding matrix (either a subsampled randomized DCT (SRTT), a Gaussian embedding matrix or a hashed randomided DCT (HRTT)).

A numerical rank estimation scheme
The numerical rank estimation scheme that follows from the observations 1. σ i (AX) ≈ σ i (A) for an approximately low-rank matrix and leading singular values, and The algorithm requires an upper bound for the numerical rank as an input (r 1 ).This input informs the size of first embedding X, and is the total number of singular values estimates that will be obtained.If the input rank r 1 was not large enough, i.e. r 1 < rank ϵ (A), Algorithm 1 is naturally only able to detect that r 1 is a lower bound for rank ϵ (A).To get a proper estimate, we then need to rerun the algorithm with a larger input rank, taking advantage of the preceding computation by appending extra columns to X and rows to Θ.
The choice of r 1 can be based on a number of possible considerations.For example, one might choose r 1 based on memory requirement, e.g. when one is unwilling to store more than (m+n)r 1 numbers.If some information is available on A (e.g. it is updated from a matrix with known singular values), that can provide a good choice of r 1 .If no information is available one might choose r 1 = O(1), say r 1 = 64; however, this can result in some inefficiency. 2  2 If the sketches employed are Gaussian, the complexity is O(mnr 1 ) (dominated by sketching) and the overhead of underestimating r 1 is minimal, bounded by a factor 2 (assuming the update rule r 1 := 2r 1 ; other increments are of course possible).However, with other sketches such as SRTT, increasing r 1 without resketching the whole matrix requires some care: For the SRTT SF D, where S is a subsampling matrix, F is a FFT/DCT matrix, and D is a diagonal matrix of random signs, we store F DAX 1 , F DAX 2 , etc., and then subsample the big matrix F DAX.
Provided the algorithm runs sucessfully, the computational complexity is O(mn log n + r 3  1 ) for a dense matrix A ∈ F m×n , where the two terms come from steps 3 and 7, assuming both X and Y are SRTT matrices.The cost of step 3 can be lower if A has structure, in which case the O(mr 1 log m) cost in step 6 may become significant.
Regarding step 7, computing the singular values of Y AX, one can of course compute them via the SVD.We do this in most cases as the O(r 3  1 ) cost is usually negligible.Alternatively, since it suffices to retrieve estimates for σ i (Y AX), one can compute the QR factorization Y AX = QR and take the diagonal elements of R as approximations of the singular values; this works as a QR factorization of a matrix whose (right) singular vectors are randomized are rank-revealing with high probability [5].

Algorithm 1: Rank Estimation with approximate orthogonalization.
Result: Given an m × n matrix A, a tolerance ϵ (optional) and an estimate for the upper bound for rank r 1 , this scheme computes an estimate for the ϵ-rank of A. If no such r exists, repeat the algorithm with a larger r 1 , e.g.r 1 := 2r 1 , by appending to the sketches.If (an estimate of) ∥A∥ is unavailable, use σ 1 (ΘAX) to approximate ∥A∥.If no tolerance ϵ is provided, plot σ i (Y AX) to detect a gap visually or consider measures as r = max i σ i (Y AX)/σ i+1 (Y AX).

Rank estimation as part of randomized low-rank approximation
It is worth discussing the case where Algorithm 1 is used as an initial step for randomized low-rank approximation, such as RSVD [19] or generalized Nyström [36].In such cases, computing the sketch AX is needed anyway, so step 3 incurs no additional cost.
Moreover, with generalized Nyström, even Y AX is required, as the low-rank approximant is of the form AX(Y AX) † Y A. Furthermore, to evaluate (Y AX) † a QR factorization Y AX = QR is computed.It follows that the additional work needed to execute Algorithm 1 is just to estimate the singular values of R. As described above, to do so one can simply take the diagonal elements of R, requiring just O(r 1 ) work.Thus in the context of generalized Nyström, a rank estimate can be obtained essentially for free.
In the next section we explain how this rank estimate can be employed to speed up fixed-precision low-rank approximation.In particular, we note that using the singular value estimates we can reduce the size of a sketch as appropriate before performing the dominant computational work.

Fixed-precision low-rank approximation problem
Rank estimation as part of low-rank approximation is inherently linked to the fixed-precision problem, where one aims to compute a low-rank factorization Â of unspecified rank such that ∥A − Â∥ ≤ ϵ∥A∥ for some user-specified tolerance ϵ.We propose a scheme to combine rank estimation and low-rank approximation through existing error estimates.The algorithm will be used here for approximation in the Frobenius norm but can easily be extended to approximation in any other norm for which randomized low-rank approximation error estimates are available.The aim will be to find an approximation of the form Â = QB, where Q ∈ F m×r has orthonormal columns and r ≤ n.
In their landmark paper of 2011, Halko, Martinsson and Tropp [19] describe and extensively analyze the randomized rangefinder (or randomized SVD) algorithm to compute such an approximation.This algorithm is displayed as Algorithm 2. Note that it can be extended to the randomized SVD algorithm by computing the SVD of B = Û ΣV T and then taking U = Q Û .A limitation of this factorization routine, as well as of other randomized low-rank approximation algorithms, is that we often wish to solve a fixed precision problem instead of a fixed rank algorithm, where the appropriate choice of the input rank r is usually unknown.
Result: Given an m × n matrix A, a target rank r ≥ 2 and an oversampling parameter p ≥ 2 such that k + p ≤ min(m, n), this scheme computes a QB-factorization of A. 1 Draw n × (r + p) random embedding matrix X. 2 [Q, ∼] = qr(AX, 0).(thin QR factorization) As mentioned, we propose to use rank estimation to tackle this problem using a priori error estimates.Under the conditions of Algorithm 2 with a Gaussian embedding matrix, it is known that [19,Thm. 10.5] We can use the Rank Estimation Algorithm 1 to compute estimates of the trailing singular values which appear in the estimate (18).Suppose we have estimates σi = σ i (Y AX) ≈ σ i (A) for the first r 1 singular values of A, as a result of step 6 of Algorithm 1. Conservatively, and taking into account that there could be underestimation in the last few estimates if the tail is light, we set σi = σr1 for i = r 1 + 1, . . ., n.We can then choose the target rank r in the rangefinder algorithm to be the first r such that where p is fixed to be 5 or 10 as suggested in [19] or where p is fixed relative to r, say p = 0.1r, as suggested in [36].
As mentioned before, the rank estimation part of this scheme requires only O(mr 1 log m + r 3  1 ) additional cost.More importantly, the singular value estimates it provides allow us to reduce the size of the sketch AX ahead of the rangefinder process.This could result in considerable speed-up, as the QR factorization of AX and the matrix-matrix multiplication Q * A are usually the dominant costs in the algorithm.Finally, we also know the rank of QB is nearly minimal while (approximately) satisfying the accuracy criteria; by contrast, reducing the rank of a QB factorization without rank estimation requires computing the SVD of B and truncating it.Result: Given an m × n matrix A, a tolerance ϵ, an oversampling parameter p ≥ 2 (either fixed absolutely or fixed relatively to r) and an upper bound for the numerical rank r 1 , this scheme computes a QB factorization of A such that ∥A − QB∥ F ≲ ϵ∥A∥.

Randomized Rangefinder:
5 Set r to be the smallest integer s.t.
2. Polynomial decay: each singular value takes the value σ i = i −p for a parameter p.We define a matrix A SP with a slow polynomial decaying spectrum (SP) with p = 1, and a matrix A F P with fast polynomial decay (FP) with p = 3.
3. Exponential decay: the singular values are logarithmically equally spaced between 1 and σ n , that is σ i = 10 −q(i−1) for a parameter q.We study a matrix A SE with slow exponential decay (SE) for which q = 0.01 and a matrix A F E with fast exponential decay (FE) for which q = 0.5.
Each of the matrices has spectral norm ∥A∥ = 1 so the relative accuracy notation will be suppressed in this section.In each of the experiments in subsections 6.2-6.4,we will choose the embedding matrix X to be a hashed randomized DCT matrix and Y = Θ to be a subsampled randomized DCT matrix.We display the perfomance of Algorithms 1 and 3 in various numerical experiments.

Detecting gaps in the spectrum
Here we demonstrate that Algorithm 1 works extremely well for approximately low-rank matrices that have gaps in their singular value spectrum.The algorithm can in fact be used in two ways: 1) if the specified tolerance is within the gap between two singular values, the algorithm will perform very well in identifying the correct rank even for a small number of samples, and 2) by considering the gaps in the singular value spectrum of Y AX, one can identify gaps in the singular value spectrum of A and use this to inform the ϵ value and target rank.
To illustrate this we compare the singular values σ i (Y AX) to σ i (A) for various values of r 1 (see Algorithm 1).Note that we oversample by 10% and disregard the singular values associated with oversampling.The results are shown in Figure 5.The figure displays the remarkable effectiveness of the algorithm, where gaps in the spectrum are clearly identified even for small values of r 1 .
We use Gaussian embeddings for the column space and SRCTs for the row space, to allow for an adaptive procedure.That is, we start with a Gaussian X of dimension 110, compute Y AX and its singular values, and then increase the sketch dimension of X by a 100 until we can identify the gap.This is how Algorithm 1 could be used if there is no upper bound on r 1 available.We Figure 5: On the left, the singular values of AG and Y AGX, 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 r1 = 110, 210, 310, 410 respectively.The gaps that exist in the spectrum of A are also visible in the spectrum of Y AGX, even for values of r1 close to the location of the gaps.On the right, the density of state plots resulting from the Chebyshev-based [42] and Lanczos-based algorithms [43].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.
compare against the performance of the DOS algorithms [42] (Chebyshev-based) and [43] (Lanczos-based).Our rank estimation algorithm took 18.33 seconds to reveal the gap at the 400th singular value.The runtime of the DOS algorithms was 2223 sec.for [42] and 199.1 sec.for [43].The DOS plots are also shown in Figure 5.
It is inherently difficult to detect gaps visible on a log-scale using a DOS algorithm.The algorithms are based on estimating a sum of Dirac delta functions using a sum of smooth functions such as Gaussians.To be able to distinguish between very small eigenvalues, the width of the Gaussian has to be as small as the smallest eigenvalue.That is, it would be necessary to have a very highresolution spectral density.The resolution refers to the size of the interval in which we estimate the number of eigenvalues [27].However, a fine resolution implies the number of sample points needed to detect larger eigenvalues is large.As a result, the computational effort necessary is not tractable.Figure 5 shows that even if the DOS is estimated perfectly for the chosen Gaussian standard deviation, the gaps are not visible.

Robustness of the rank estimator
The other four example matrices discussed (SP, FP, SE, and FE) do not have significant gaps in their singular value spectra, which makes the numerical rank less well-defined.Rather than focusing on the numerical rank estimate, r, it thus makes sense to look at σ r+1 (A) and its proximity to ϵ. Figure 6 shows the results of Algorithm 1 applied 100 times to each of the example matrices.We see that, although the rank estimate r we recover is often an underestimate when there are no significant gaps between singular values near the desired accuracy ϵ, σ r+1 (A) is always very close to ϵ.One could argue this is not necessarily a shortcoming of the algorithm, but rather a display of the difficulties associated to the concept of the exact numerical rank rank ϵ (A) when σ r (A) ≈ σ r+1 (A).Importantly, the algorithm always found a rank that satisfies the two desiderata set out in Section 1.1.In particular, when there are significant gaps in the singular value spectrum we recover the exact numerical rank in all except one instance; see the rightmost column of the figure (and the forthcoming Figure 5).

Comparison with other numerical rank estimation algorithms
Suppose that the goal is to obtain the numerical rank of a given matrix for some tolerance ϵ.We compare five different methods to achieve this for a

Discussion
Let us conclude with a few extensions and open problems related to Algorithm 1.
As presented, Algorithm 1 does not take into account the structure of A such as symmetry or sparsity.It would be of interest to devise a variant that respects such structure.
When the input rank r 1 is insufficient, i.e. r 1 < rank ϵ (A), Algorithm 1 needs to be rerun with a larger r 1 .It would be helpful (if possible) to be able to find an appropriate new r 1 from the information obtained from the first run.
While we primarily focused on rank estimation, the analysis indicates that many of the (leading) singular values of A can be estimated by σ i (ΘAX).Tuning each σ i (ΘAX) to obtain a reliable (ideally unbiased) estimator is left for future work.

Figure 1 :
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 σ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 σr+1 = σ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 :
Figure2: One could argue that in none of the cases discussed in Figure1is the numerical rank very well-defined, because there is no clear gap between any of the singular values.We repeat the experiment of Figure1but 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 :
Figure 3: Numerical rank estimation of matrices with polynomially decaying singular values and exponentially decaying singular values, where only the sketch AX ∈ R m×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, rankϵ(A).

Figure 4 :
Figure 4: Relative error in the singular value estimates of a real matrix B with fast polynomially decaying singular values of size 10 5 × 2000 based on the singular values of Y B, where Y ∈ R 4000×10 5is an embedding matrix (either a subsampled randomized DCT (SRTT), a Gaussian embedding matrix or a hashed randomided DCT (HRTT)).

Algorithm 3 :
Rank Estimation combined with Randomized Rangefinder.

Figure 6 :
Figure6: The results of Algorithm 1 applied 100 times to four dense synthetic test matrices with singular value spectra that decay Slow Polynomially (SP), Fast Polynomially (FP), Slow Exponentially (SE), and Fast Exponentially (FE), as defined at the start of Section 6.Each matrix is square diagonal with n = 10 5 .The blue line indicates the tolerance ϵ and the blacks dots represent the singular values σi(A).The red dots indicate the singular values σ r+1 (A), where r is the numerical rank estimate found by the algorithm.The area of the dot represents the frequency with which it a specific r was found.The algorithm would work perfectly if the singular value on the blue line were found in each iteration, as is (almost) the case for the FE matrix.The upper and lower row correspond to upper limits r1 = 2 rankϵ(A) and r1 = 4 rankϵ(A) respectively.

Figure 8 :
Figure 8: We compare four algorithms to compute a low-rank QB approximation to solve the fixed precision problem ∥ASE − QB∥ ≤ ϵ, where ASE is a 10 5 × 10 5 square matrix with slow exponential decaying singular values.The algorithms include rank estimation with rangefinder (RE), rangefinder (RF) and RandQB_EI.We compare for the parameter choices r1 = 1500, 3000, r = 1500, 3000 and b = 100 for each of the algorithms respectively.

Figure 9 :
Figure9: We compare three algorithms for fixed precision estimation that require only two views of the matrix, randomized rangefinder (Algorithm 2), randomized SVD with truncation and rank estimation combined with randomized rangefinder (Algorithm 3, shown as RE (QB)).For a fixed tolerance ϵ = 10 −3 , we vary the input target rank/upper bound for the rank.Although each of the algorithms satisfy the condition of an error smaller than ϵ, we see that Algorithm 3 is much faster, especially for greater overestimates.The matrix is ASE of dimension n = 10 5 .