low-rank approximation, matrix-free algorithm, adaptive rank determination, adjoint-free range finder, fixed-precision problem, randomized norm estimation

1 Introduction↩︎

We consider the fixed-tolerance low-rank approximation problem. Given a matrix \(\mathbf{A} \in \mathbb{R}^{m \times n}\) and a target tolerance \(\varepsilon > 0\), we seek to compute matrices \(\mathbf{Q} \in \mathbb{R}^{m \times k}\) and \(\mathbf{B} \in \mathbb{R}^{k \times n}\) such that \[\label{eq:auv} \| \mathbf{A} - \mathbf{Q} \mathbf{B} \| \le \tau,\tag{1}\] where the rank \(k \ll \min(m,n)\) is a priori unknown and must be determined dynamically. In 1 , the tolerance \(\tau\) is a user-prescribed parameter, specified either as an absolute threshold \(\tau = \varepsilon\) or a relative threshold \(\tau = \varepsilon \, \|\mathbf{A}\|\). We focus on the Frobenius norm \(\|\cdot\|_\text{F}\) and the spectral norm \(\|\cdot\|_2\). As is typical of randomized methods, the algorithms we develop satisfy 1 with high probability and, in the spectral norm, up to a modest constant (4).

We operate in the matrix-free setting: \(\mathbf{A}\) is accessible only through matrix-vector products (matvecs)—the forward and adjoint actions \(\mathbf{x} \mapsto \mathbf{A}\mathbf{x}\) and \(\mathbf{y} \mapsto \mathbf{A}^\top\mathbf{y}\)—while its individual entries are not. This setting is ubiquitous in computational science. In Bayesian inverse problems and uncertainty quantification, for instance, the Hessian of the log-posterior is defined implicitly, yet a low-rank approximation of it accelerates downstream optimization and sampling [1], [2]. There, a single matvec may require a full forward and adjoint PDE solve, so the dominant cost is the number of operator applications, not the internal dense linear algebra. Our goal is therefore to solve 1 while minimizing the number of matvecs, and to do so accurately for \(\varepsilon\) down to machine precision.

We focus on a specific structural form of 1 , the QB approximation [3], [4], where \(\mathbf{Q}\) has orthonormal columns; the choice \(\mathbf{B} = \mathbf{Q}^\top\mathbf{A}\) then minimizes the residual in the Frobenius norm. QB is a building block for many matrix factorizations, decoupling range identification (the range finder) from entry approximation (the projection): examples include the truncated SVD [5], the interpolative decomposition [6], the higher-order SVD for tensor compression [7], and hierarchical matrix approximations [8], [9]. Each relies on the fixed-tolerance QB primitive that the matrix-free methods of 3 and 4 provide.

1.1 Existing Work and Limitations↩︎

1.1.0.1 Classical deterministic methods

Direct factorization methods such as the singular value decomposition (SVD) and rank-revealing QR (RRQR) [10], [11] offer excellent stability and deterministic error bounds, but are designed for explicitly stored matrices: they require entry-wise access and scale as \(\mathcal{O}(\min(m,n)mn)\), so they are inapplicable in the matrix-free setting.

Iterative Krylov subspace methods such as Lanczos and Arnoldi [12] use only matvecs to extract the leading singular triplets. Although effective for rapidly decaying spectra, they are numerically unstable in finite precision—the Lanczos vectors lose orthogonality, requiring expensive re-orthogonalization [13]—and converge slowly for flat or clustered spectra, requiring many sequential matvecs.

Cross and skeleton methods such as the adaptive cross approximation (ACA) [14], [15] build skeleton decompositions from a sublinear number of rows and columns. Although scalable for smooth kernels (e.g., boundary element methods), ACA must evaluate individual entries \(A_{ij}\) at arbitrary positions; for implicit operators such as Hessians, a single entry costs a full matvec, negating this advantage. Its greedy pivoting is also prone to instability or premature convergence on matrices with zero blocks, noise, or localized features [16].

1.1.0.2 Modern randomized methods

Randomized numerical linear algebra now provides a mature toolkit for low-rank approximation [17], [18]. Randomized rank-revealing factorizations, such as randomized column-pivoted QR [19], [20], use random embeddings to select pivot columns efficiently, bypassing the communication overhead of classical pivoting. However, they still update the trailing matrix explicitly, so they cannot operate in a matrix-free setting.

Fixed-rank matrix-free methods such as the randomized SVD (RSVD) [21], the Nyström method for positive semi-definite matrices [22], and the generalized Nyström scheme [23] are matrix-free and parallelizable, but fixed-rank: the user must supply the target rank a priori (which can itself be estimated by randomized numerical-rank estimation [24]). Pass-efficient variants further reduce the number of views (passes) over \(\mathbf{A}\) [25], but they too require the rank in advance. In multiscale applications such as \(\mathcal{H}\)-matrix construction [8], [26]–[28], where the numerical rank fluctuates across thousands of off-diagonal blocks, a fixed rank causes either memory over-allocation or loss of fidelity.

To address adaptivity, Halko, Martinsson, and Tropp [21] introduced an adaptive range finder, hereafter randQB_HMT, the state-of-the-art (SOTA) adaptive method in the spectral norm; it is based on a probabilistic error bound that is pessimistic in practice and over-estimates the rank. Yu, Gu, and Li [4] introduced randQB_EI, the SOTA adaptive method in the Frobenius norm, which builds the range block by block and tracks the residual error dynamically. However, it estimates the residual Frobenius norm by Gramian-type tracking (\(\mathbf{A}\mathbf{A}^\top\)), squaring the effective condition number and capping the attainable precision near \(\sqrt{\varepsilon_{\text{mach}}}\) (\(\approx 10^{-8}\) in IEEE double precision). Because its indicator subtracts \(\|\mathbf{B}\|_{\mathrm{F}}^2\) from \(\|\mathbf{A}\|_{\mathrm{F}}^2\), it also requires the Frobenius norm \(\|\mathbf{A}\|_{\mathrm{F}}\) to be computed explicitly from \(\mathbf{A}\) and is thus not matrix-free.

A related adaptive Frobenius-norm method, by Gorman et al. [29], attains comparable accuracy but certifies the error by post-processing the random sketches from all iterations, making it substantially more expensive than the rank pruning introduced here (3.2). More recently, Liu and Yu [30] proposed adaptive, fixed-threshold randomized algorithms with error analysis in both the Frobenius and spectral norms; like the methods above, however, they operate on the explicit entries of \(\mathbf{A}\) and are not matrix-free. Adaptive randomized schemes have likewise been developed for other factorizations and settings, including rank-adaptive CUR [31], adaptive low-rank approximation for image compression [32], and dynamical low-rank approximation [33].

In short, no existing method is simultaneously matrix-free, adaptive in rank, and accurate down to machine precision in both the Frobenius and spectral norms; closing this gap is the goal of the present work.

1.2 Contributions↩︎

The primary contribution is a family of adaptive, matrix-free QB algorithms that solve the fixed-tolerance low-rank approximation problem down to machine precision (\(\varepsilon_{\text{mach}}\)) in both the Frobenius and spectral norms. They build on the blocked randomized QB iteration [3], [4], which constructs \(\mathbf{Q}\) and \(\mathbf{B}\) incrementally in blocks of \(b\) columns per iteration. Our main developments are:

• Randomized residual indicators: Rather than a single, generic indicator, we develop a family of randomized error indicators tailored to each method—the sketched residual norm in the Frobenius case and a look-ahead block norm in the spectral case, each with an adjoint-free counterpart. Every indicator is read off the block sketch already computed at each iteration, at negligible cost, and avoids the failure modes of prior adaptive schemes: the “precision wall” near \(\sqrt{\varepsilon_{\text{mach}}}\) of energy-subtraction (Gramian) indicators [4], and the rank over-estimation caused by the conservative probabilistic bound of [21], remaining accurate and tight down to machine precision.

• Rank Pruning and the Block-Size Dilemma: Our most significant algorithmic innovation is a post-iteration rank pruning step. Large block sizes \(b\) are needed for matrix-matrix (BLAS-3) efficiency, but they over-estimate the rank when the tolerance is met mid-block. Our pruning decouples the block size from the final rank, preserving high performance while retaining a near-minimal rank.

• Adjoint-Free Adaptivity: We propose adjoint-free variants that determine the rank without the adjoint \(\mathbf{A}^\top\). When the adjoint is expensive or unavailable and only the basis \(\mathbf{Q}\) is needed, these variants roughly halve the passes over the data (\(\mathbf{B}\) is not formed).

We emphasize that the basic matrix-free and adjoint-free blocked QB frameworks have appeared in prior work; our contributions are the rank-pruning step, the matrix-free estimate of \(\|\mathbf{A}\|\) for relative tolerances, and a thorough numerical study of the resulting methods.

1 summarizes the four resulting methods, organized by norm and by whether the adjoint operator is required.

2 Preliminaries and Problem Formulation↩︎

2.1 Notation↩︎

Throughout this work, matrices are denoted by bold uppercase letters (e.g., \(\mathbf{A}, \mathbf{B}\)), vectors by bold lowercase letters (e.g., \(\mathbf{x}, \mathbf{y}\)), and scalars by standard lowercase or Greek letters (e.g., \(k, \varepsilon\)). The singular value decomposition (SVD) of \(\mathbf{A}\) is denoted as \(\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^\top\), with singular values \(\sigma_1 \ge \sigma_2 \ge \dots \ge 0\). We utilize two standard matrix norms to measure the approximation error:

• The Frobenius norm, defined as \(\|\mathbf{A}\|_{\text{F}} = \sqrt{\sum_{i,j} |A_{ij}|^2} = \sqrt{\sum_j \sigma_j^2}\).

• The spectral norm, defined as \(\|\mathbf{A}\|_2 = \sigma_1\).

Our pseudocode uses two MATLAB-style primitives: \(\texttt{randn}(n, b)\) returns an \(n \times b\) matrix with independent standard Gaussian \(\mathcal{N}(0,1)\) entries, and \(\texttt{orth}(\mathbf{Y})\) returns a matrix whose columns form an orthonormal basis for the range of \(\mathbf{Y} \in \mathbb{R}^{n \times b}\) with \(b < n\), computed in practice via a thin (unpivoted) QR factorization.

2.2 The Standard Randomized QB Decomposition↩︎

The standard randomized algorithm for computing a low-rank QB approximation proceeds in two primary phases: range finding and projection. The method comes with a priori error guarantees: the approximation error is within a polynomial factor of the optimal value \(\sigma_{k+1}(\mathbf{A})\) with high probability [21]; see [34] for a unified error analysis in both the Frobenius and spectral norms. It does, however, require the target rank \(k\) to be fixed in advance. The randomized QB method goes as follows:

1. Range Finding: Draw a Gaussian random test matrix \(\mathbf{\Omega} \in \mathbb{R}^{n \times (k+p)}\) (Here, \(p\) is an oversampling parameter chosen typically as 5 or 10). Compute the sketch \(\mathbf{Y} = \mathbf{A}\mathbf{\Omega}\) and form an orthonormal basis for its range: \(\mathbf{Q} = \texttt{orth}(\mathbf{Y})\).

2. Projection: Project the original matrix onto this basis to form the core matrix: \(\mathbf{B} = \mathbf{Q}^\top\mathbf{A}\).

For adaptive, fixed-tolerance settings, the blocked QB method [3], [4] builds \(\mathbf{Q}\) and \(\mathbf{B}\) incrementally in blocks of size \(b\). At the \(i\)-th iteration, let \(\mathbf{Q}^{(i-1)}\) and \(\mathbf{B}^{(i-1)}\) be the accumulated bases. The algorithm sketches the residual operator \(\mathbf{E}_{i-1} = \mathbf{A} - \mathbf{Q}^{(i-1)}\mathbf{B}^{(i-1)}\) as \(\mathbf{Y}_i = \mathbf{E}_{i-1} \mathbf{\Omega}_i\) using a block of random vectors \(\mathbf{\Omega}_i \in \mathbb{R}^{n \times b}\). Once \(\mathbf{Q}_i = \texttt{orth}(\mathbf{Y}_i)\) is computed, the block \(\mathbf{B}_i\) is formed via \(\mathbf{B}_i = \mathbf{Q}_i^\top\mathbf{A}\). This process is summarized in 1

Notice that 1 measures the error at each iteration via the (Frobenius or spectral) norm of the residual matrix \(\mathbf{E}_i\). Forming \(\mathbf{E}_i\) explicitly violates the matrix-free constraint: it costs as many matvecs as \(\mathbf{A}\) has rows or columns.

2.3 Estimating the Frobenius Norm↩︎

We recall a classical result on randomized Frobenius-norm estimation, due to Hutchinson [35], which underlies our residual tracking without forming matrix entries.

1 establishes unbiasedness in expectation; concentration bounds [36] further guarantee that a single realization stays close to the true norm with high probability. For an order-of-magnitude estimate of the residual error, a small block size (\(b = 5\) or \(10\)) suffices, so convergence can be monitored with cheap, low-dimensional sketches.

2.4 Estimating the Spectral Norm↩︎

For stopping criteria in the spectral norm, we must track the singular values of the residual (and of \(\mathbf{A}\)) from sketched information alone. We review a bound relating the singular values of the projection \(\mathbf{B}= \mathbf{Q}^\top\mathbf{A}\) to those of \(\mathbf{A}\) in the QB framework of 2.2.

By 2, the spectral norm of the projection \(\mathbf{B}\) is a reliable proxy for the dominant singular values of \(\mathbf{A}\): with \(i=1\), \(\|\mathbf{A}\|_2 = \sigma_1(\mathbf{A}) \le \eta_1 \sigma_1(\mathbf{B}) = \eta_1 \|\mathbf{B}\|_2\). Applying this with the residual \(\mathbf{E}_i = \mathbf{A}- \mathbf{Q}^{(i)}\mathbf{B}^{(i)}\) in place of \(\mathbf{A}\) lets us monitor the spectral-norm error without ever forming the residual.

A second, complementary tool is the notion of a subspace embedding [17]: a sketch that preserves norms on a subspace also preserves its singular values. The simplest such embedding is a (scaled) random Gaussian matrix [17], which is the sketch used throughout this work.

3 Adaptive Matrix-free Algorithms in Frobenius norm↩︎

Suppose \(\mathbf{A}\in \mathbb{R}^{m \times n}\) is accessible only through its forward and adjoint actions \(\mathbf{x}\mapsto \mathbf{A}\mathbf{x}\) and \(\mathbf{y}\mapsto \mathbf{A}^\top \mathbf{y}\). Our objective is to compute \(\mathbf{Q}\in \mathbb{R}^{m \times k}\) with orthonormal columns and \(\mathbf{B}= \mathbf{Q}^\top\mathbf{A}\in \mathbb{R}^{k \times n}\) such that \[\|\mathbf{A}- \mathbf{Q}\mathbf{B}\|_{\mathrm{F}} \leq \tau,\] where the target rank \(k\) is not known in advance and must be determined dynamically. The choice \(\mathbf{B}= \mathbf{Q}^\top\mathbf{A}\) minimizes the Frobenius-norm error for a given basis \(\mathbf{Q}\). The tolerance is specified either as an absolute threshold \(\tau = \varepsilon\) or a relative threshold \(\tau = \varepsilon \, \|\mathbf{A}\|_{\mathrm{F}}\), for some \(\varepsilon > 0\).

Standard approaches, such as the blocked QB algorithm [3], explicitly form the residual matrix \(\mathbf{E}_i\) at every iteration (see the residual update in 1). Consequently, they cannot be applied in the matrix-free setting. In this section, we demonstrate how to adapt the blocked QB framework into a purely matrix-free algorithm by resolving the following key algorithmic challenges:

1. Residual Estimation (3.1): We introduce a randomized error indicator that efficiently estimates the residual norm \(\|\mathbf{E}_i\|_{\mathrm{F}}\) without forming \(\mathbf{E}_i\). This bypasses the need for explicit matrix entries, enabling an adaptive stopping criterion relying exclusively on matrix-vector products.

2. Block Size Selection (3.2): We analyze strategies for choosing the block size parameter to maximize computational efficiency while preventing severe overestimation of the numerical rank.

3. Adjoint-Free Variant (3.3): We introduce a forward-only modification that eliminates the need for adjoint queries (\(\mathbf{y}\mapsto \mathbf{A}^\top \mathbf{y}\)) when only the basis \(\mathbf{Q}\) is required.

3.1 Randomized Error Indicator and Matrix-Free Algorithm↩︎

Recall the challenge in the blocked QB method is that the stopping criterion requires the error \(\|\mathbf{E}_{i}\|_{\mathrm{F}}\), where the residual matrix \[\label{eq:ei} \mathbf{E}_{i} = \mathbf{A}- \sum_{k=1}^{i}\mathbf{Q}_k\mathbf{B}_k \quad i=1,2,\ldots.\tag{2}\] To avoid accessing entries of \(\mathbf{E}_{i}\), we use Hutchinson’s trace (Frobenius norm) estimator. But recall that a sample matrix \(\mathbf{Y}_i = \mathbf{E}_{i-1} \, \mathbf{\Omega}_i\) is already computed to build the basis at every iteration in the blocked QB algorithm. Here, \(\mathbf{E}_{i-1}\) is the residual from the previous steps, and \(\mathbf{\Omega}_i\) is a newly drawn random test matrix. Throughout, we normalize the random blocks so that \(\mathbf{\Omega}_i \in \mathbb{R}^{n \times b}\) has i.i.d. \(\mathcal{N}(0, 1/b)\) entries—equivalently, \(\mathbf{\Omega}_i = \texttt{randn}(n,b)/\sqrt{b}\)—which absorbs the \(1/\sqrt{b}\) factor that would otherwise appear in every norm estimate. Applying 1 to \(\sqrt{b}\,\mathbf{\Omega}_i\) then dictates that \[\begin{align} \mathbb{E} \left[ \|\mathbf{Y}_i\|_{\mathrm{F}}^2 \right] =\mathbb{E} \left[ \|\mathbf{E}_{i-1}\mathbf{\Omega}_i\|_{\mathrm{F}}^2 \right] =\|\mathbf{E}_{i-1}\|_{\mathrm{F}}^2. \end{align}\] Therefore, we propose using \(\|\mathbf{Y}_i\|_{\mathrm{F}}^2\) as an unbiased randomized estimator of the current residual norm (squared) \(\|\mathbf{E}_{i-1}\|_{\mathrm{F}}^2\) (see 1).

Importantly, this estimator requires only the product \(\mathbf{E}_{i-1} \mathbf{\Omega}_i\) and does not require explicit access to the entries of the residual. Furthermore, evaluating this estimator incurs negligible \(\mathcal{O}(bm)\) additional computational cost.

We no longer need \(\mathbf{E}_{i-1}\) for evaluating the error, but we still need to compute \(\mathbf{Y}_i = \mathbf{E}_{i-1} \, \mathbf{\Omega}_i\). How do we do this without storing \(\mathbf{E}_{i-1}\)? Assuming \(i-1\) iterations are complete, we use a simple associativity trick: \[\begin{align} \label{Y95i32-32full32description} \mathbf{Y}_i =\mathbf{E}_{i-1}\mathbf{\Omega}_i & =\left(\mathbf{A}-\sum_{k=1}^{i-1}\mathbf{Q}_k\mathbf{B}_k\right)\mathbf{\Omega}_i \\ & = \underbrace{\mathbf{A}\mathbf{\Omega}_i}_{\text{black-box Mat-Vec}} -\sum_{k=1}^{i-1}\mathbf{Q}_k \underbrace{(\mathbf{B}_k\mathbf{\Omega}_i)}_{\text{small matrix mult}}. \notag \end{align}\tag{3}\] Forming \(\mathbf{Y}_i\) this way completely eliminates the need to ever form or store the residual matrix \(\mathbf{E}_{i-1}\).

The combination of the randomized error indicator and the implicit residual expansion gives rise to a fully adaptive, matrix-free algorithm, summarized in 2. A defining feature of this new framework is that the sample matrix \(\mathbf{Y}_i\) fulfills a dual algorithmic role. First, in its traditional capacity, it serves as the subspace generator to construct the next orthonormal block \(\mathbf{Q}_i\). Second, by simply evaluating its Frobenius norm \(\|\mathbf{Y}_i\|_{\mathrm{F}}\), it simultaneously provides an estimate of the current residual error.

We now highlight two practical considerations for the robust and efficient implementation of 2.

1. Loss of Orthogonality: In exact arithmetic, the output matrix \(\mathbf{Q}\) naturally maintains orthonormal columns. In practice, finite-precision arithmetic inevitably leads to a gradual loss of orthogonality. To counteract this, a re-orthogonalization step against the previously computed basis blocks is employed after evaluating the new basis vectors.

2. Memory Allocation: The iterative algorithm computes the output matrices \(\mathbf{Q}\) and \(\mathbf{B}\) one block at a time. Dynamically reallocating memory for these matrices with increasing sizes at every iteration incurs a severe performance penalty. An optimized implementation should instead pre-allocate memory up to a maximum expected rank. If the iterative rank exceeds this estimate, the array sizes can be doubled to amortize the cost of reallocation.

3.2 Block Size Selection and Rank Pruning↩︎

Modern cache-based computing architectures are most efficient when computations can be organized as matrix-matrix products of sufficient size. For this reason, choosing a sufficiently large block size \(b\) in 2 can substantially improve practical performance. The drawback is that the output rank grows in multiples of \(b\), meaning a large block size may overestimate the true numerical rank by as much as \(b-1\), leading to unnecessary storage and downstream computational costs.

To mitigate this, we introduce a lightweight rank pruning procedure that prunes the final block produced by 2. The pruning operates exclusively on the last computed block (the last \(b\) rows of \(\mathbf{B}\)), and it has negligible computational overhead, just \(\mathcal{O}(b n)\) operations. The objective is to keep the computed low-rank approximation as compact as possible while still rigorously satisfying the prescribed tolerance.

Suppose 2 is applied to \(\mathbf{A}\in\mathbb{R}^{m\times n}\) and the look-ahead test passes at iteration \(i \ge 2\), so the algorithm retains \(\tilde{\mathbf{Q}} = [\mathbf{Q}_1\; \mathbf{Q}_2 \; \cdots\; \mathbf{Q}_{i-1} ]\) and \(\tilde{\mathbf{B}} = [\mathbf{B}_1^\top\; \mathbf{B}_2^\top \; \cdots\; \mathbf{B}_{i-1}^\top ]^\top\). The state at termination reads: \[\begin{align} \text{Previous step failed:} \quad &{\|\mathbf{E}_{i-2}\|_\mathrm{F} > \tau} \quad \text{where} \;\; \mathbf{E}_{i-2} = \mathbf{A}- \sum_{k=1}^{i-2}\mathbf{Q}_k \mathbf{B}_k \\ \text{Current step passed:} \quad &{\|\mathbf{E}_{i-1}\|_\mathrm{F} \le \tau} \quad \;\; \text{where} \;\; \mathbf{E}_{i-1} = \mathbf{E}_{i-2} - \mathbf{Q}_{i-1}\mathbf{B}_{i-1}. \end{align}\] Here, \(\mathbf{Q}_{i-1}\in\mathbb{R}^{m\times b}\) and \(\mathbf{B}_{i-1}\in\mathbb{R}^{b\times n}\) are the final blocks added to the basis. (Technically speaking, we don’t have the exact errors \(\|\mathbf{E}_{i-2}\|_\mathrm{F}\) and \(\|\mathbf{E}_{i-1}\|_\mathrm{F}\) but their randomized estimators.)

To formalize this approach, we partition the index set \(\{1,\ldots,b\}\) into two disjoint subsets: \[\{1,\ldots,b\} = \mathbf{s}\cup\mathbf{r}, \qquad \mathbf{s}\cap\mathbf{r}=\emptyset, \qquad |\mathbf{s}|+|\mathbf{r}|=b.\] Here, \(\mathbf{s}\) denotes the indices to select and \(\mathbf{r}\) denotes the indices to remove. Let \(\mathbf{Q}_{i-1}^{\mathbf{s}}\) and \(\mathbf{Q}_{i-1}^{\mathbf{r}}\) contain the corresponding columns in \(\mathbf{Q}_{i-1}\), and let \(\mathbf{B}_{i-1}^{\mathbf{s}}\) and \(\mathbf{B}_{i-1}^{\mathbf{r}}\) contain the corresponding rows in \(\mathbf{B}_{i-1}\). It follows that the final block update can be split as \(\mathbf{Q}_{i-1}\mathbf{B}_{i-1} = \mathbf{Q}_{i-1}^{\mathbf{s}}\mathbf{B}_{i-1}^{\mathbf{s}} + \mathbf{Q}_{i-1}^{\mathbf{r}}\mathbf{B}_{i-1}^{\mathbf{r}}.\) Selecting the smallest necessary set \(\mathbf{s}\) corresponds to solving \[\begin{align} \label{1st32opt32prob} \min\;|\mathbf{s}| \quad\text{subject to}\quad \big\| \mathbf{E}_{i-2} - \mathbf{Q}_{i-1}^{\mathbf{s}}\mathbf{B}_{i-1}^{\mathbf{s}} \big\|_{\mathrm{F}} \le \tau. \end{align}\tag{5}\] Once \(\mathbf{s}\) is determined, the compact form of our low-rank approximation becomes \({\mathbf{Q}} = [\mathbf{Q}_1\; \cdots\; \mathbf{Q}_{i-2} \; \mathbf{Q}_{i-1}^{\mathbf{s}} ]\) and \({\mathbf{B}} = [\mathbf{B}_1^\top\; \cdots\; \mathbf{B}_{i-2}^\top \; (\mathbf{B}_{i-1}^{\mathbf{s}})^\top ]^\top\), which satisfy \(\|\mathbf{A}- {\mathbf{Q}} {\mathbf{B}}\|_{\mathrm{F}} \le \tau\) with high probability.

Solving 5 is equivalent to solving the following “dual” formulation: \[\begin{align} \label{2nd32opt32prob} \max\;|\mathbf{r}| \quad\text{subject to}\quad \big\| \mathbf{E}_{i-1} + \mathbf{Q}_{i-1}^{\mathbf{r}}\mathbf{B}_{i-1}^{\mathbf{r}} \big\|_{\mathrm{F}} \le \tau. \end{align}\tag{6}\] This equivalence follows immediately from the fact that \(|\mathbf{r}|=b-|\mathbf{s}|\) and \[\begin{align} \mathbf{E}_{i-1} + \mathbf{Q}_{i-1}^{\mathbf{r}}\mathbf{B}_{i-1}^{\mathbf{r}} &= \mathbf{E}_{i-1} + \left( \mathbf{Q}_{i-1}\mathbf{B}_{i-1} - \mathbf{Q}_{i-1}^{\mathbf{s}}\mathbf{B}_{i-1}^{\mathbf{s}} \right) \\ &= \left( \mathbf{E}_{i-1} + \mathbf{Q}_{i-1}\mathbf{B}_{i-1} \right) - \mathbf{Q}_{i-1}^{\mathbf{s}}\mathbf{B}_{i-1}^{\mathbf{s}} \\ &= \mathbf{E}_{i-2} - \mathbf{Q}_{i-1}^{\mathbf{s}}\mathbf{B}_{i-1}^{\mathbf{s}} . \end{align}\]

Computational Challenge: A naive approach for solving 6 is to start with \(\|\mathbf{E}_{i-1}\|_\mathrm{F}\) and increase it by sequentially adding rank-\(1\) outer products. However, forming the intermediate matrices explicitly is computationally expensive, and our method would lose its matrix-free advantage.

Instead, we only need to track how the Frobenius norm changes when we add rank-\(1\) outer products from the final block \(\mathbf{Q}_{i-1} \mathbf{B}_{i-1}\). Because the columns of \(\mathbf{Q}_{i-1}\) are orthonormal, adding or subtracting a rank-1 component only changes the overall Frobenius norm by the magnitude of the row vector in \(\mathbf{B}_{i-1}\). This is stated as the following theorem:

Consequently, applying 3 to the final block, the optimization problem in 6 reduces to selecting the largest subset of rows from \(\mathbf{B}_{i-1}\) whose squared Frobenius norms can be added to the squared residual norm without violating \(\|\mathbf{E}_{i-1}\|_{\mathrm{F}}^{2} + \|\mathbf{B}_{i-1}^{\mathbf{r}}\|_{\mathrm{F}}^{2} \le \tau^{2}.\) To summarize, we do not need to form any explicit matrices.

To retain as few rank-1 outer products as possible, we sort the row norms of matrix \(\mathbf{B}_{i-1}\) in ascending order; start with \(\|\mathbf{E}_{i-1}\|_\mathrm{F}^2\); and increase it by sequentially adding the row norms (squared) to just cross \(\tau^2\). This is summarized in 3. Of course, we do not have the exact \(\|\mathbf{E}_{i-1}\|_\mathrm{F}\) (and \(\|\mathbf{A}\|_\mathrm{F}\)) but instead use their randomized estimators in [eq:est95me] 4 .

3.3 Adjoint-free Adaptive Range Finder↩︎

For applications where only the orthonormal basis \(\mathbf{Q}\) is required and the projection \(\mathbf{B}= \mathbf{Q}^\top\mathbf{A}\) is not explicitly needed, we present a variant of 2 that requires only the forward operator \(\mathbf{A}\). By eliminating the need for the adjoint operator \(\mathbf{A}^\top\), this variant is both more broadly applicable and computationally faster. This task is commonly referred to as the range finder problem [17], [21].

In 2, the matrix \(\mathbf{B}\) is utilized in two specific locations: (i) the computation of the implicit residual \(\mathbf{Y}_i\) via the term \(\mathbf{B}_k \mathbf{\Omega}_i\), and (ii) the rank-pruning step which requires the row norms of \(\mathbf{B}_{i-1}\). We address these two requirements without the adjoint operator as follows.

First, we observe that the term \(\mathbf{B}_k \mathbf{\Omega}_i\) (\(k=1,2,\ldots, i-1\)) can be evaluated by shifting the parentheses to the right: \[\label{eq:bi} \mathbf{B}_k \mathbf{\Omega}_i = (\mathbf{Q}_k^\top \mathbf{A}) \mathbf{\Omega}_i = \mathbf{Q}_k^\top (\mathbf{A}\mathbf{\Omega}_i) = \mathbf{Q}_k^\top \mathbf{Z}_i\tag{7}\] By computing \(\mathbf{Q}_k^\top \mathbf{Z}_i\), we obtain the necessary projection using only the already-computed sample matrix \(\mathbf{Z}_i\).

Second, the rank-pruning procedure typically requires the row norms of the final block \(\mathbf{B}_{i-1}\) to determine which columns of \(\mathbf{Q}_{i-1}\) are redundant. In the adjoint-free case, we propose using the row norms of the sketch \(\tilde{\mathbf{B}}_{i-1} = \mathbf{B}_{i-1} \mathbf{\Omega}_i\) as a proxy. According to 1, the squared row norms of \(\tilde{\mathbf{B}}\) are unbiased Hutchinson estimators of the squared row norms of \(\mathbf{B}_{i-1}\). Specifically, for the \(j\)-th row of the final block \(\mathbf{B}_{i-1}\), denoted by \(\mathbf{b}_j\), the pruning algorithm relies on the estimator \(\hat{w}_j = \|\mathbf{b}_j \mathbf{\Omega}_i\|_2^2\). While we have established that \(\mathbb{E}[\hat{w}_j] = \|\mathbf{b}_j\|_2^2\), the reliability of the pruning decision depends on the variance of this estimator, which is given by \(\text{Var}(\hat{w}_j) = \frac{2}{b} \|\mathbf{b}_j\|_2^4\). This variance decreases as the block size \(b\) increases, meaning that larger block sizes not only improve BLAS-3 efficiency but also enhance the precision of the adjoint-free rank refinement. In practice, even with moderate block sizes (e.g., \(b=16\)), the probability of significantly misidentifying the “low-energy” rows is small. Furthermore, since the pruning step is only applied to the final block, any slight over-retention of columns due to estimator variance results only in a tiny increase in the final rank.

We summarize our adjoint-free method in 4.

4 Adaptive Matrix-Free Algorithms in Spectral Norm↩︎

We now extend the matrix-free framework to the spectral norm (operator 2-norm). In this setting, the objective is to determine a target rank \(k\) dynamically such that the following approximation holds: \[\label{eq:aqb95spec} \| \mathbf{A} - \mathbf{Q}\mathbf{B} \|_2 \le \tau,\tag{8}\] where \(\mathbf{Q} \in \mathbb{R}^{m \times k}\) has orthonormal columns and \(\mathbf{B} = \mathbf{Q}^\top\mathbf{A} \in \mathbb{R}^{k \times n}\). The tolerance \(\tau\) is a user-defined threshold, provided as either an absolute value \(\tau = \varepsilon\) or a relative value \(\tau = \varepsilon \|\mathbf{A}\|_2\).

As in the Frobenius-norm case, we seek to construct \(\mathbf{Q}\) and \(\mathbf{B}\) incrementally using only matrix-vector products with \(\mathbf{A}\) and \(\mathbf{A}^\top\). However, the spectral norm presents a unique challenge: How to estimate the spectral norm of the residual matrix given a pair of factors \(\mathbf{Q}\) and \(\mathbf{B}\). To address this, we resolve the following key algorithmic challenges:

1. Residual Estimation (4.1): We introduce a “look-ahead” randomized error indicator. By utilizing the fact that the QB framework naturally targets the dominant singular spectrum, the norm of the most recently computed block \(\mathbf{B}_i\) provides an efficient estimate for the residual norm of the preceding approximation.

2. Block Size Selection (4.2): We analyze strategies for choosing the block size parameter to maximize computational efficiency while preventing severe overestimation of the numerical rank.

3. Adjoint-Free Variant (4.3): We introduce a forward-only modification that eliminates the need for adjoint queries (\(\mathbf{y}\mapsto \mathbf{A}^\top \mathbf{y}\)) when only the basis \(\mathbf{Q}\) is required.

4.1 Randomized Error Indicator and Matrix-Free Algorithm↩︎

The primary difficulty in an adaptive spectral-norm QB algorithm is that the stopping criterion depends on \(\|\mathbf{E}\|_2\), where \(\mathbf{E} = \mathbf{A} - \mathbf{Q} \mathbf{B}\) for a given pair of factors \(\mathbf{Q}\) and \(\mathbf{B}\). We observe that at each iteration, the blocked QB method, i.e., 1, generates a basis \(\mathbf{Q}_i\) whose span approximates the subspace spanned by dominant left singular vectors of the residual \(\mathbf{E}_{i-1}\). If the block size \(b\) is sufficiently large (e.g., \(b=16\)), then, with probability almost 1, the error \(\| \mathbf{E}_{i-1} - \mathbf{Q}_i \mathbf{Q}_i^\top \mathbf{E}_{i-1} \|_2\) is close to the minimum error in rank-\((b-p)\) approximation [21], where \(p\) is the oversampling parameter (say \(p=5\) or \(p=10\)).

Consequently, the matrix \[\mathbf{B}_i = \mathbf{Q}_i^\top \mathbf{A} = \mathbf{Q}_i^\top \left( \mathbf{E}_{i-1} + \sum_{k=1}^{i-1} \mathbf{Q}_k \mathbf{B}_k \right) = \mathbf{Q}_i^\top \mathbf{E}_{i-1}\] captures the action of the residual on its most significant subspace, where the first equality is by the definition of \(\mathbf{B}_i\), the second equality is by the definition of \(\mathbf{E}_{i-1}\) in 2 , and the last equality holds by the orthogonality of the basis matrices \(\mathbf{Q}_i\).

We summarize our adaptive, matrix-free algorithm in 5.

4.2 Block Size Selection and Rank Pruning↩︎

The block-size dilemma described in 3.2 applies equally here: a large \(b\) improves BLAS-3 throughput but lets the output rank \(k = i \cdot b\) overshoot the minimal rank needed to meet the tolerance. We again prune the final retained block, but now in the spectral norm: instead of discarding rows by their norm, we truncate \(\mathbf{B}_{i-1}\) according to its singular value distribution.

The change of criterion is dictated by the norm. In the Frobenius case, the orthogonality of \(\mathbf{Q}_{i-1}\) makes the squared error additive over the rows of \(\mathbf{B}_{i-1}\) (3), so pruning reduces to dropping rows of small norm. The spectral norm is not additive in this way; the quantity that governs it is the singular-value spectrum of the block \(\mathbf{Q}_{i-1} \mathbf{B}_{i-1}\). Since \(\mathbf{Q}_{i-1}\) has orthonormal columns, we compute the SVD of \(\mathbf{B}_{i-1}\) and, invoking the Eckart–Young theorem, retain exactly the leading \(r\) singular directions whose singular values reach the tolerance, discarding the trailing directions whose singular values—and hence whose contribution to the residual—fall below \(\varepsilon\).

The procedure for rank pruning is summarized in 6. Since \(\mathbf{B}_{i-1}\) has dimensions \(b \times n\), its SVD costs only \(\mathcal{O}(nb^2)\) flops—negligible relative to the matvecs—yet it keeps the output rank close to the optimal numerical rank of the operator. The following theorem makes the tolerance guarantee precise.

4.3 Adjoint-free Adaptive Range Finder↩︎

As in the Frobenius-norm case (3.3), many applications require only the orthonormal basis \(\mathbf{Q}\) and not the projection \(\mathbf{B}= \mathbf{Q}^\top\mathbf{A}\). We therefore present a variant of 5 that accesses \(\mathbf{A}\) through the forward operator alone, eliminating the adjoint queries \(\mathbf{y}\mapsto \mathbf{A}^\top\mathbf{y}\).

In 5 the adjoint is used in two places: (i) the look-ahead indicator \(\|\mathbf{B}_i\|_2\), where \(\mathbf{B}_i = \mathbf{Q}_i^\top\mathbf{A}\), and (ii) the rank-pruning step, which truncates \(\mathbf{B}_{i-1}\) by its singular values. We replace both with forward-only quantities. Both replacements rest on the same principle that 1 justifies: a matrix’s singular values are read from a Gaussian sketch of it.

First, recall that \(\mathbf{E}_{i-1} = \mathbf{A}- \sum_{k=1}^{i-1}\mathbf{Q}_k\mathbf{Q}_k^\top\mathbf{A}\) denotes the residual after the first \(i-1\) blocks. Its random sketch \(\mathbf{E}_{i-1}\mathbf{\Omega}_i\) can be assembled using only the forward action \(\mathbf{Z}_i = \mathbf{A}\mathbf{\Omega}_i\): \[\label{eq:af95resid95spec} \mathbf{E}_{i-1}\mathbf{\Omega}_i = \left(\mathbf{A}- \sum_{k=1}^{i-1}\mathbf{Q}_k\mathbf{Q}_k^\top\mathbf{A}\right)\mathbf{\Omega}_i = \mathbf{Z}_i - \sum_{k=1}^{i-1}\mathbf{Q}_k\left(\mathbf{Q}_k^\top\mathbf{Z}_i\right) = \mathbf{Y}_i,\tag{10}\] where we have used \((\mathbf{Q}_k^\top\mathbf{A})\mathbf{\Omega}_i = \mathbf{Q}_k^\top\mathbf{Z}_i\) to avoid forming \(\mathbf{B}_k\). As a result, the spectral norm of this sketch yields a look-ahead estimate of the residual, \[\label{eq:af95indicator95spec} \|\mathbf{Y}_i\|_2 = \| \mathbf{E}_{i-1}\mathbf{\Omega}_i \|_2 \approx \|\mathbf{E}_{i-1}\|_2 ,\tag{11}\] the spectral-norm analogue of the Frobenius indicator of 3.3. It replaces \(\|\mathbf{B}_i\|_2\) as a crude stopping criterion and is estimated by a few steps of power iteration applied directly to \(\mathbf{Y}_i\).

Second, the pruning step of 6 truncates \(\mathbf{B}_{i-1}\) by its singular values, but \(\mathbf{B}_{i-1} = \mathbf{Q}_{i-1}^\top\mathbf{A}\) is never formed here. Exactly as in the Frobenius-norm case 7 , we replace it with the forward sketch \[\label{eq:af95prune95spec} \tilde{\mathbf{B}}_{i-1} = \mathbf{Q}_{i-1}^\top\mathbf{Z}_i = \mathbf{Q}_{i-1}^\top\mathbf{A}\,\mathbf{\Omega}_i = \mathbf{B}_{i-1}\mathbf{\Omega}_i ,\tag{12}\] obtained from the current block’s sample \(\mathbf{Z}_i = \mathbf{A}\mathbf{\Omega}_i\) at no additional matvec cost. Since \(\mathbf{\Omega}_i^\top\) is, with high probability, a subspace embedding for the row space of \(\mathbf{B}_{i-1}\), 1 gives \(\sigma_j(\tilde{\mathbf{B}}_{i-1}) \approx \sigma_j(\mathbf{B}_{i-1})\). Computing the SVD \(\tilde{\mathbf{B}}_{i-1} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\top\) and retaining the \(r\) left singular vectors with \(\sigma_j(\tilde{\mathbf{B}}_{i-1}) \geq \varepsilon\) gives the pruned basis \(\mathbf{Q}= [\mathbf{Q}_1\;\ldots\;\mathbf{Q}_{i-2}\; \left(\mathbf{Q}_{i-1}\mathbf{U}_r \right)]\), the adjoint-free counterpart of 6.

We summarize the procedure in 7. A caveat is in order: guaranteeing a genuine \(\varepsilon\)-subspace embedding would require the block size \(b\) to grow like \(1/\varepsilon^2\), far more than is practical, so the sketch delivers only a same-order estimate of the singular values rather than a highly accurate one. We emphasize, then, that 8 is to be understood in an approximate sense: unlike the Frobenius-norm methods of 3, which certify the prescribed bound (with high probability), the spectral-norm algorithms do not guarantee 8 exactly, but instead keep the computed error \(\|\mathbf{A} - \mathbf{Q}\mathbf{B}\|_2\) on the order of the tolerance \(\tau\).

5 Numerical Results↩︎

We evaluate the proposed matrix-free adaptive algorithms with two objectives: to show (1) that the adaptive stopping criteria accurately capture the target tolerance across a range of singular-value decay profiles, and (2) that the rank-pruning strategies yield a near-optimal output rank. All computations are performed in double precision.

5.0.0.1 Reproducibility

The MATLAB implementations of the proposed algorithms, together with the driver scripts that reproduce all numerical experiments and figures reported in this section, are publicly available at https://github.com/SMART-Algebra/Adaptive-Matfree-Lowrank.

5.1 Experimental Setup and Baselines↩︎

We compare the proposed algorithms against the following baselines:

1. Frobenius-norm baseline (randQB_EI [4]): the SOTA adaptive QB method. Its error estimate relies on a “top-down” subtraction of squared Frobenius norms, which is prone to catastrophic cancellation in finite precision and typically limits its reliable tolerance to about \(\mathcal{O}(\sqrt{\varepsilon_{\text{mach}}})\) ([rmk:randQB95ei]). We show that our matrix-free variants remove this limitation, remaining stable down to machine precision.

2. Spectral-norm baseline (randQB_HMT [21]): the SOTA adaptive randomized range finder of Halko, Martinsson, and Tropp. It uses only the forward operator and returns only the orthonormal basis \(\mathbf{Q}\), so it is adjoint-free. Its a posteriori error estimate is provably reliable but pessimistic ([21]), so it typically overestimates the rank.

3. Theoretical optimum (truncated SVD, t-SVD): for any tolerance \(\varepsilon\), the Eckart–Young theorem identifies the truncated SVD as the rank-\(k\) approximation that minimizes the error in both the Frobenius and spectral norms. The SVD cannot be applied in the matrix-free setting without first forming the entire matrix, but we include it as a “gold standard” that measures the rank over-determination factor of our algorithms — how many extra dimensions the proposed randomized methods need to match the error of the optimal (but infeasible) SVD.

5.1.0.1 Test matrices

We test our methods on matrices exhibiting a range of decay profiles, from rapid exponential decay (typical of smoothing operators) to slow algebraic decay (common in certain statistical datasets) to S-shape decay. Let \(\mathbf{A}= \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\top\), where \(\mathbf{U}, \mathbf{V}\in \mathbb{R}^{1024 \times 1024}\) are random Gaussian matrices and \(\mathbf{\Sigma} \in \mathbb{R}^{1024 \times 1024}\) is diagonal with entries defined below:

• Matrix 1: \(\sigma_i = \text{exp}(-i/2)\) \(1\leq i \leq 1024\) Exponential decay

• Matrix 2: \(\sigma_i = 1/i^p\) \(1\leq i \leq 1024\) Polynomial decay

• Matrix 3: \(\sigma_i =\begin{cases} 10^{-6(i-1)/399}, & 1 \le i \le 400, \\[6pt] 10^{\,a\, (i-400)^{-\alpha} + b}, & 401 \leq i \leq 1024, \end{cases}\) S-shape decay
  where \(\alpha = 10^{-0.8} \approx 0.158\), and the constants \(a\) and \(b\) are fixed by requiring the tail to join the first piece continuously at \(\sigma_{401} = 10^{-6}\) and to reach \(\sigma_{1024} = 5\varepsilon_{\text{mach}}\); this gives \(a \approx 14\) and \(b \approx -20\).

The two exponents for Matrix 2 are chosen to stress-test the baselines, not to constrain our methods. For the Frobenius-norm experiments we take \(p=4\), so that part of the spectrum falls below \(\sqrt{\varepsilon_{\text{mach}}}\) and thereby exposes the precision wall of randQB_EI. For the spectral-norm experiments we take \(p=2\), so that the numerical ranks are relatively large and all methods require many iterations to converge, which sharpens the contrast between the proposed methods and randQB_HMT. These choices are designed to highlight the differences between the methods and are not limitations of the proposed algorithms.

In all cases, the tolerance \(\tau\) is specified as a relative threshold: \(\tau = \varepsilon \|\mathbf{A}\|,\) for a wide range of \(\varepsilon\), where \(\|\cdot\|\) denotes either the Frobenius or spectral norm, depending on the algorithm under test. The relative criterion is more demanding than an absolute one, as the algorithm must additionally estimate \(\|\mathbf{A}\|\) on the fly; our methods do so matrix-free from the same sketches ([rmk:relative32tolerance]), whereas randQB_EI computes \(\|\mathbf{A}\|_{\mathrm{F}}\) explicitly from \(\mathbf{A}\) and is thus not matrix-free.

We use a block size of \(b = 32\) throughout, except for Matrix 1, whose small numerical rank calls for \(b = 16\) to keep the iteration counts comparable across matrices. For each method we report the achieved (exact) relative error in the relevant norm. For randQB_EI and the matrix-free variants (randQB_MF_Fro, randQB_MF_Spec), it is \(\|\mathbf{A}- \mathbf{Q}\mathbf{B}\| / \|\mathbf{A}\|\), where \(\mathbf{Q}\) and \(\mathbf{B}\) are the computed pair. The adjoint-free variants (randQB_AF_Fro, randQB_AF_Spec) and randQB_HMT never form the projection \(\mathbf{B}= \mathbf{Q}^\top\mathbf{A}\), so their error is evaluated as \(\|\mathbf{A}- \mathbf{Q}\mathbf{Q}^\top\mathbf{A}\| / \|\mathbf{A}\|\). For the t-SVD it is the optimal error of the best rank-\(k\) approximation.

5.2 Comparison in the Frobenius Norm↩︎

We compare four methods: the truncated SVD (t-SVD; the optimal, but not matrix-free, baseline), randQB_EI (SOTA, but not matrix-free), and the two proposed methods randQB_MF_Fro and randQB_AF_Fro. [fig:fro_exp,fig:fro_poly,fig:fro_sshape] report the results for the exponential, polynomial, and S-shaped spectra (Matrices 1–3), respectively. Each figure shows the computed rank versus the prescribed relative tolerance \(\varepsilon\) (left) and the per-iteration error indicator of each method (right), together with a table of the computed rank \(k\) and the achieved (exact) relative error in the Frobenius norm (evaluated as described in 5.1; tabulated at every other tolerance for brevity, while the panels use the full range).

5.2.0.1 Accuracy and tolerance range

Across all three spectra, the proposed randQB_MF_Fro matches randQB_EI wherever the latter is reliable, while remaining accurate over a far wider range of tolerances. Its achieved relative error tracks the prescribed tolerance almost perfectly, staying just at or below \(\varepsilon\) all the way down to machine precision. The adjoint-free variant randQB_AF_Fro is slightly less stable but still attains an achieved error of the same order as \(\varepsilon\) throughout. In contrast, randQB_EI is reliable only for \(\varepsilon \gtrsim \mathcal{O}(\sqrt{\varepsilon_{\text{mach}}})\); below that wall it breaks down (entries marked fail in the tables).

5.2.0.2 Rank efficiency

The ranks produced by both proposed methods stay close to the optimal ranks of the t-SVD (typically within a handful of extra dimensions) so the matrix-free, randomized rank determination costs little. Unlike the t-SVD, however, the proposed methods access \(\mathbf{A}\) only through matrix-vector products and never form a dense factorization, which makes them far cheaper for large operators.

5.2.0.3 Error indicators and the breakdown of randQB_EI

The methods differ in how they monitor the residual error \(\|\mathbf{A} - \mathbf{Q}\mathbf{B}\|_{\mathrm{F}}\). randQB_EI uses the identity \(\|\mathbf{E}\|_{\mathrm{F}}^2 = \|\mathbf{A}\|_{\mathrm{F}}^2 - \|\mathbf{B}\|_{\mathrm{F}}^2\). While exact in infinite precision, this “top-down” subtraction is highly susceptible to catastrophic cancellation once the residual is small relative to \(\|\mathbf{A}\|_{\mathrm{F}}\), which (as noted in [4]) restricts the reliable tolerance to roughly \(\varepsilon \geq 2.7 \times 10^{-7}\) in double precision. Below this wall the method breaks down in one of two ways: either (i) the tracked quantity \(\|\mathbf{A}\|_{\mathrm{F}}^2 - \|\mathbf{B}\|_{\mathrm{F}}^2\) turns negative, so the indicator is undefined; or (ii) the indicator stops decreasing near \(\mathcal{O}(\sqrt{\varepsilon_{\text{mach}}})\) and the stopping criterion is never met, so the iteration runs to its maximum number of steps and returns a full-rank factorization with no low-rank approximation. Which behavior occurs is driven by the random sketch rather than the spectrum: either can arise for any test matrix, and the same matrix may exhibit one or the other across runs; the figures report one representative failure each. In contrast, randQB_MF_Fro (and randQB_AF_Fro) use the randomized indicator \(\|\mathbf{Y}\|_{\mathrm{F}} = \| ( \mathbf{A} - \mathbf{Q}\mathbf{B}) \mathbf{\Omega} \|_{\mathrm{F}} \approx \|\mathbf{A} - \mathbf{Q}\mathbf{B}\|_{\mathrm{F}},\) where \(\mathbf{\Omega}\) is a normalized Gaussian matrix (i.i.d.\(\mathcal{N}(0,1/b)\) entries) with \(b\) columns. Being a norm of a sketched residual, this indicator is inherently non-negative and free of the cancellation that plagues the energy-subtraction approach, so it reliably drives the iteration down to tolerances near machine precision.

5.3 Comparison in the Spectral Norm↩︎

We now extend the comparison to the spectral norm, evaluating four methods: the t-SVD (the optimal, but not matrix-free, baseline), the SOTA range finder randQB_HMT [21], and the two proposed methods randQB_MF_Spec and randQB_AF_Spec. [fig:spec_exp,fig:spec_poly,fig:spec_sshape] report the results for the exponential, polynomial, and S-shaped spectra (Matrices 1–3), respectively. Each figure shows the computed rank versus the prescribed relative tolerance \(\varepsilon\) (left) and the per-iteration error indicator of each method at the smallest tolerance (right), together with a table of the computed rank \(k\) and the achieved (exact) relative error in the spectral norm (evaluated as described in 5.1; tabulated at every other tolerance for brevity, while the panels use the full range).

5.3.0.1 Rank overestimation and comparison with randQB_HMT

A significant advantage of our method is its rank efficiency compared to the existing adaptive range finder. The existing method utilizes an error indicator that includes a safety constant \(C = 10\sqrt{2/\pi}\). While this constant provides a rigorous probabilistic bound, it acts as a significant “amplification” factor that causes the algorithm to perceive the error as much larger than it actually is in many cases.

As illustrated in the left panels of 11 12 13, this pessimistic stopping criterion leads randQB_HMT to consistently overestimate the rank required to satisfy the tolerance — often inflating it to (near-)full rank at the tightest tolerances. In contrast, randQB_MF_Spec and randQB_AF_Spec avoid this over-determination, yielding numerical ranks that are far more compact — close to the optimal t-SVD ranks for the rapidly-decaying spectra, and well below randQB_HMT’s throughout.

5.3.0.2 Accuracy of the spectral error indicator

A key feature of randQB_MF_Spec is the use of the look-ahead block norm \(\|\mathbf{B}_i\|_2\) as a proxy for the total residual error. As shown in the right panels of 11 12 13, this indicator provides a remarkably tight estimate of the true residual error throughout the iteration process. The indicator remains reliable across many orders of magnitude.

6 Conclusions↩︎

We have introduced a family of adaptive, matrix-free randomized QB algorithms for high-precision low-rank approximation. Together, they address three bottlenecks of existing methods: excessive data passes, rank over-determination, and the precision wall inherent in Gramian-based error estimation. Our numerical experiments demonstrate several advantages of the proposed framework:

• Robustness in the Frobenius Norm: While the existing method randQB_EI fails to satisfy tolerances below \(\mathcal{O}(\sqrt{\varepsilon_{\text{mach}}})\), our matrix-free error indicator reliably tracks the residual error down to machine precision.

• Efficiency in the Spectral Norm: By utilizing a look-ahead block norm indicator, we avoid the pessimistic constants found in an existing adaptive range finder, yielding significantly more compact ranks that closely match the optimal truncated SVD.

• Reduced Data Movement: When only the orthonormal basis is required, the adjoint-free variants dispense with the backward pass over the data, roughly halving the number of operator applications while degrading the approximation accuracy only slightly.

Several extensions are natural:

6.0.0.1 Tensor compression

Inserting the relative-tolerance range finder into the sequentially truncated HOSVD [7], [40] yields a matrix-free ST-HOSVD for the Tucker format. Standard randomized implementations estimate the per-mode truncation error with a Gramian-based indicator of the randQB_EI [4] type, which is limited to tolerances of order \(\mathcal{O}(\sqrt{\varepsilon_{\text{mach}}})\) ([rmk:randQB95ei]); our indicators track the residual to machine precision, so the resulting ST-HOSVD overcomes this precision wall.

6.0.0.2 Hierarchical matrices

The construction of hierarchical (\(\mathcal{H}\)-)matrices, currently under active development, is a closely related target. Randomized black-box constructions build the representation from matrix–vector products, either by peeling off the off-diagonal blocks level by level [9], [41], [42] or, more recently, by simultaneously compressing and factorizing the matrix [43]. In all of these, the numerical rank fluctuates unpredictably across the many off-diagonal blocks, so fixed-rank range finders either over-allocate memory or lose fidelity, whereas the adaptive, fixed-tolerance methods of 3 and 4 set the rank block-by-block to a prescribed accuracy.

6.0.0.3 High-performance software

Because the algorithms are inherently blocked, they map naturally onto BLAS-3 kernels and modern parallel hardware, making a high-performance (GPU and distributed-memory) implementation a natural next step, building on recent HPC randomized-factorization efforts [20], [44]–[47]. A further direction is to extend these matrix-free techniques to streaming data.

Acknowledgments↩︎

The authors used AI-based writing assistants to help edit language and improve the readability of the manuscript. All technical content was conceived, verified, and approved by the authors, who take full responsibility for the work.

References↩︎

---

1. Department of Mathematics, North Carolina State University, Raleigh, NC ().↩︎

2. Department of Mathematics, North Carolina State University, Raleigh, NC ().↩︎

3. Department of Mathematics, North Carolina State University, Raleigh, NC ().↩︎