1 Introduction↩︎

Given a matrix \(A\in\mathbb{R}^{m\times n}\) a common task is to compute \(k\) “good” columns of \(A\). Let \(\mathcal{C}\) denote a subset of \([n]\) with \(\lvert \mathcal{C}\rvert = k\). Two simple criteria we may seek to optimize are the error approximating \(A\) by a rank-k matrix constrained to use \(A(:,\mathcal{C})\) as a basis for the column space, i.e., \[\label{eq:lrapprox} \|A-P_\mathcal{C}A\|_{\{2,F\}},\tag{1}\] where \(P_\mathcal{C}\) is the orthogonal projector onto the span of \(A(:,\mathcal{C})\), and the quality of the basis, i.e., \[\label{eq:smin} \sigma_{\min}{(A(:,\mathcal{C}))}.\tag{2}\] The latter criterion is, e.g., particularly relevant when searching for basic solutions to least-squares problems. It is in this context that Golub, Klema, and Stewart introduced their scheme for selecting \(\mathcal{C}\) [1]; henceforth we refer to their scheme as GKS.

GKS proceeds in two stages. First, the leading \(k\) right singular vectors of \(A\) are computed and we denote the associated matrix \(V_k\in\mathbb{R}^{n\times k}\). Then, a column-pivoted QR factorization [2], [3] is applied to \(V_k^T\) to select \(\mathcal{C}\). The first bounds for this method in terms of the quality of the pivoted QR factorization can be found in [4] for the two norm and [5] for the Frobenius norm. The key point is that GKS reduces the column based low-rank approximation and basis conditioning problems to computing a rank-revealing QR factorization of a matrix with orthonormal rows.

When computing the necessary pivoted QR factorization of \(V_k^T\) one has two choices. The common choice is an algorithm due to Golub and Businger [2], [3]. Implementations of this algorithm are widely available and the performance is often good in practice. However, its worst case theoretical performance is abysmal (with bounds that grow exponentially in \(k\)). In fact, recent work [6] provides pathological examples (for sufficiently large matrices) that exhibit poor behavior even under the orthonormal row constraint—addressing an open problem [7]. If good theoretical guarantees are desired, algorithms that compute strong rank-revealing QR factorizations (such as that of Gu and Eisenstat [8]) can be used. However, these algorithms can be much less efficient in practice (particularly as the desired bounds are strengthened).

A pivoting scheme due to Stewart [9] (based on work by Bischof [10]) is preferable to both options in this context. We refer to the scheme as Bischof-Stewart pivoting. In 3 we show that when applied to a matrix with orthonormal rows Bischof-Stewart pivoting computes a strong rank-revealing QR factorization with stronger theoretical guarantees than provided by the Gu and Eisenstat algorithm when limited to polynomial computational complexity. Then, in 5 we provide a randomized variant of the pivoting strategy for matrices with orthonormal rows that provides the same theoretical guarantees while being much faster in practice. While the runtime is not deterministic, it strongly outperforms (e.g., by up to two orders of magnitude) the deterministic variant in all cases we consider.

After the introduction of necessary background material in 2 and the analysis of Bischof-Stewart pivoting in 3, 4 studies the effectiveness of GKS when augmented with Bischof-Stewart pivoting, which we refer to as modified GKS (mGKS). In the \(\lVert \cdot \rVert_2\) norm, we show that for 1 mGKSis as good as a strong rank-revealing QR factorization (see, e.g., [8]); the same is true for 2 . Specifically, we show that mGKSselects \(\mathcal{C}\) such that \[\lVert A-P_\mathcal{C}A \rVert_2 \leq \sqrt{1+k(n-k)}\sigma_{k+1}(A)\] and \[\sigma_{\min}(A(:,\mathcal{C}))\geq \frac{\sigma_k(A)}{\sqrt{1+k(n-k)}}.\] When \(m\geq n\) this algorithm provides the same theoretical guarantees as recent work by Osinsky [11] in \(\lVert \cdot \rVert_2\) while being much more computationally efficient since excessive residual computations are avoided.² In the \(\lVert \cdot \rVert_F\) mGKSmay (theoretically) lag behind Osinsky’s algorithm for 1 , but using results from [5] we are able to characterize matrices where the bounds are similar (though, again, mGKSis cheaper).

To address situations where \(V_k\) is not available and is, instead, approximated by some \(\hat{V}_k,\) 4 highlights how in either the two or Frobenius norm mGKSsimply pays a sub-optimality factor of \[\frac{\lVert A-A\hat{V}_k\hat{V}_k^T \rVert_{\{2,F\}}}{\lVert A-AV_k V_k^T \rVert_{\{2,F\}}},\] which naturally captures how much worse \(\hat{V}_k\) is at approximating \(A\) than \(V_k\). Notably, \(\lVert A-A\hat{V}_k\hat{V}_k^T \rVert_{\{2,F\}}\) provides a suitable target when there is no appreciable gap between \(\sigma_k\) and \(\sigma_{k+1}\) and approximation of the dominant subspace is ill-posed [12]. This means (approximate) mGKSis well suited to situations where \(A\) is only available via matrix-vector products.³ In this setting there are many ways to compute \(\hat{V}_k\) but there is no natural way to compute a strong rank-revealing QR factorization of \(A\) directly.

Column subset selection is simultaneously a classical problem in numerical linear algebra, and one that merits, and receives, continued attention. 6 highlights several ways in which the results herein can be extended (e.g., to two sided interpolatory factorizations), improved (e.g., to try and close the theoretical gap with volume-sampling), and provide natural improvements in application areas.

1.1 Related work↩︎

In the Frobenius norm, Cortinovis and Kressner [13] provide a randomized algorithm that, in expectation, achieves the same upper bound as Osinsky without residual computations. Leveraging the fact that the algorithm in [13] is equivalent to \(k\)-volume sampling of a projection matrix, Epperly [14] gives an independent proof of the same approximation bounds and provides more efficient sampling methods. The origins of these methods are in the (k-)volume sampling literature [15]–[19], which we will not survey here.

This manuscript provides an improved understanding of how to leverage randomness within column selection algorithms. At one extreme sampling strategies can be highly efficient (e.g., uniform sampling), albeit with weak error guarantees. Conversely, computing a (full) strong rank-revealing factorization or performing “exact” k-volume sampling is quite expensive. A natural middle ground is to either approximate \(V_k\) or “sketch” \(A\) as \(GA\) using some appropriate \(G\) and then apply a column selection algorithm to \(V_k^T\) or \(GA\) (e.g., a simple pivoted QR). Recent work [20], [21] elaborates on this tradeoff.

At least in the two norm, the most natural blend of theoretical guarantees and practical performance is to leverage randomness to approximate \(V_k\) and then follow that up with (randomized) Bischof-Stewart pivoting applied to \(V_k^T\). To get strong theoretical guarantees from \(GA\) typically requires appealing to strong rank-revealing factorizations [22], [23], which could be inefficient in practice. Random sampling schemes with strong error guarantees such as those highlighted earlier [13], [14] provide no better guarantees than mGKSin the two norm (regardless of whether true or approximate singular vectors are used). Moreover, their error guarantees are probabilistic (i.e., in expectation)—those of mGKSare always satisfied, even when our randomized pivoting scheme is used.

2 Background↩︎

2.1 Matrix Factorizations↩︎

We will primarily leverage two key matrix factorizations: (partial) pivoted QR factorizations and the singular value decomposition (SVD). Given a permutation matrix \(\Pi\in\mathbb{R}^{n\times n}\) and any \(k\leq\min{(m,n)}\) we let \[\label{eq:pivotedQR} A\Pi = \begin{bmatrix}Q_1 & Q_2\end{bmatrix}\begin{bmatrix}R_{11} & R_{12} \\ & M\end{bmatrix}\tag{3}\] denote a partially computed (full) QR factorization of \(A\Pi\) where \(Q_1\in\mathbb{R}^{m\times k}\) and \(Q_2\in\mathbb{R}^{m\times (m-k)}\) have orthonormal columns, \(R_{11}\in\mathbb{R}^{k\times k}\) is upper triangular, \(M\in\mathbb{R}^{(m-k) \times (n-k)}\), and \(R_{12}\in\mathbb{R}^{k\times (n-k)}\). Observe that, since we place no structural restrictions on \(M,\) we have defined a partial reduction to upper triangular form for \(k<\min{(m,n)}\); for this reason we use \(M\) rather than \(R_{22}\). Lastly, it will sometimes be convenient to split the permutation \(\Pi\) as \[\Pi = \begin{bmatrix}\Pi_1 & \Pi_2\end{bmatrix},\] where \(\Pi_1\) is the first \(k\) columns of \(\Pi\) and \(\Pi_2\) the remaining \(n-k\). When \(k=\min{(m,n)}\) we call this a pivoted QR factorization of \(A\).

For the SVD, again given any \(k\leq\min{(m,n)}\) we denote \[\label{eq:svd} A = U_k\Sigma_k V_k^T + U_\perp \Sigma_\perp V_\perp^T,\tag{4}\] where \(U_k\in\mathbb{R}^{m\times k},\) \(V_k\in\mathbb{R}^{n\times k},\) and \(\Sigma_k\in\mathbb{R}^{k\times k}\) are used to denote the leading \(k\) left singular vectors, right singular vectors, and singular values, and \(U_\perp\in\mathbb{R}^{m\times (m-k)},\) \(V_\perp\in\mathbb{R}^{n\times (n-k)},\) \(\Sigma_\perp\in\mathbb{R}^{(m-k)\times (n-k)}\) denote the remaining singular values and vectors. As usual, \(U_k,V_k,U_\perp,\) and \(V_\perp\) have orthonormal columns and \(\Sigma_k\) and \(\Sigma_\perp\) are diagonal with non-negative entries denoted \(\sigma_1,\ldots,\sigma_{\min{(m,n)}}\).

2.2 (Strong) rank-revealing QR factorizations↩︎

It follows from 3 that \(\lVert A-P_\mathcal{C}A \rVert = \lVert M \rVert\) for any unitarily invariant norm and \(\sigma_{\min}(A(:,\mathcal{C})) = \sigma_{\min}(R_{11})\). Therefore, the questions we asked in the introduction are intrinsically tied to pivoted QR factorizations. Accordingly, there has been significant work on algorithms that compute \(\Pi\) to ensure that the QR factorization 3 has certain desirable properties. Of particular note is the development of (strong) rank-revealing QR factorizations (see, e.g., [4], [24] for a discussion of rank-revealing and the refinement to strong rank-revealing in [8]).

For a factorization to be strong rank-revealing it must satisfy stringent criteria. Rather than just requesting that \(\|M\|_2\approx \sigma_{k+1}\) and \(\sigma_{\min}(R_{11}) \approx \sigma_k,\) the singular values of \(R_{11}\) must well approximate \(\sigma_1,\ldots,\sigma_k,\) the singular values of \(M\) must well approximate \(\sigma_{k+1},\ldots,\sigma_{\min{(m,n)}},\) and the interpolation coefficients \((R_{11}^{-1}R_{12})_{ij}\) must have controlled size. Precise statements may be found in [8]; recently, [25] provided necessary and sufficient conditions for \(\Pi\) to induce a strong rank-revealing factorization.⁴

This manuscript is primarily concerned with the reduced criteria in 1 and 2 . Nevertheless, we will show that the modified pivoting strategy we discuss next computes a strong rank-revealing factorization in certain settings. Moreover, strong rank-revealing factorizations give good bounds for those problems and provide a natural point of comparison. From a computational standpoint, a key consideration is that algorithms for strong rank-revealing QR typically have to allow for “swaps.” In other words, for some fixed \(k\) a candidate \(\Pi_1\) is computed (e.g., via the standard pivoted QR algorithm of Golub and Businger [2], [3]) and then “fixed” by computing swaps between \(\Pi_1\) and \(\Pi_2\) till the desired bounds are achieved. Such a scheme is less efficient (albeit marginally so if the number of swaps is small) and a key consideration is how to bound the number of swaps—doing so necessitates slightly relaxing the most stringent criteria that could be asked for. For mGKSwe need no such relaxation.

2.3 Stewart’s modified pivoting criterion↩︎

7 contains a comprehensive description of the Bischof-Stewart pivoting scheme along with numerical experiments to illustrate practical performance. Here we briefly summarize the relevant ideas. Given a partially computed QR factorization as in 3 the core idea is to pick the next pivot (i.e., column \(k+1\) of \(\Pi\)) by choosing the \(j\) that minimizes \[\label{eq:BSPivotBrief} \frac{1+\|R_{11}^{-1}R_{12}e_j\|_2^2}{\|Me_j\|_2^2}.\tag{5}\] By convention, at the first step (i.e., \(k=0\) to \(k=1\)) we simply choose the \(j\) that maximizes \(\|Me_j\|_2^2\) (i.e., the largest column of \(A\)). The standard Golub-Businger pivoting strategy simply sticks to this selection criterion at each step, always choosing \(j\) to maximize \(\|Me_j\|_2^2\).

The motivation behind minimizing 5 is to minimize the growth in \(\|R_{11}^{-1}\|_F\). To see this, consider a partially computed QR factorization of \(A\in\mathbb{R}^{m\times n}\) after \(k \leq \min{(m,n)}\) steps of the form 3 . If we let \(\hat{R}_{11}\) denote the upper left \((k+1) \times (k+1)\) block of 3 after column \(j\) of the block matrix \[\begin{bmatrix} R_{12} \\ M\end{bmatrix}\] is chosen as a pivot (i.e., columns \(k+1\) and \(k+j\) of \(A\Pi\) are swapped) and the factorization is progressed by one column (i.e., one additional Householder reflector is applied after the swap to reduce column \(k+1\) of \(A\Pi\) to upper triangular form), then \[\label{eq:RinvF95update} \|\hat{R}_{11}^{-1}\|_F^2 = \|R_{11}^{-1}\|_F^2 + \frac{1+\|R_{11}^{-1}R_{12}e_j\|_2^2}{\|Me_j\|_2^2}.\tag{6}\] This is a simple computation that follows from the fact that \[\hat{R}_{11}^{-1} = \begin{bmatrix} R_{11}^{-1} & -\frac{R_{11}^{-1}R_{12}e_j}{\|Me_j\|_2} \\ & \frac{1}{\|Me_j\|_2} \end{bmatrix}.\]

This small change has, in certain circumstances, profound impact on the theoretical performance. 7.2 shows that such a scheme can be nearly as computationally efficient as the Golub-Businger scheme.⁵ Moreover, its form enables a fast randomized variant with strong theoretical guarantees. A natural question, therefore, is why the modified pivoting strategy didn’t “take off.” While impossible to provide a definitive answer, one possible reason is that for certain pathological problems such as the Kahan matrix [26] the Bischof and Stewart scheme also fails. Its value is only theoretically evident when applied to matrices with orthonormal rows—as in the GKS framework.

3 Bischof-Stewart is strong rank revealing for matrices with orthonormal rows↩︎

Given a matrix with orthonormal rows 1 characterizes how well the pivoting scheme of Bischof and Stewart controls singular values of the selected subset. Notably, Osinsky [11] proved the \(i=m\) case in the analysis of that paper’s main algorithmic contribution (which is distinct from mGKS), and Avron and Boutsidis [27] provide similar guarantees albeit with an algorithm that is quadratic in the number of columns (though that cost can be reduced [28]). 1 highlights the connection to Osinsky’s work and that this algorithm achieves the quality of a maximum-volume submatrix [29], albeit without any guarantees of having maximal volume. For 1 we provide a direct proof based on 6 , though aspects of the strategy mirror Osinsky’s work.

1 shows that when applied to matrices with orthonormal rows Bischof-Stewart pivoting computes a strong rank-revealing QR factorization (per [8] and [25]) for \(i = 1,\ldots,m\). To formalize this statement we extract from the proof of 1 that at step \(i\) we have \(\|R_{11}^{-1}R_{12}\|_{\max} \leq \sqrt{1 + \frac{i(n-m)}{m-i+1}},\) and observe that when \(i<m\) \(\sigma_{\ell}(M) = \sigma_{\ell+i}(Z^T)\).

In fact, in the setting of orthonormal rows 1 shows that Bischof-Stewart pivoting is preferable to the strong rank-revealing algorithm of Gu and Eisenstat [8]. Bischof-Stewart yields stronger control on \(\sigma_i(R_{11})\) for \(i<m,\) and to get as strong of a result Algorithm 4 of [8] at \(i=m\) would require setting \(f=1\) (at which point the algorithm’s runtime may not be polynomial in the size of the matrix).

4 Theoretical guarantees for mGKS↩︎

We first give a brief description of the mGKSframework and then provide a theoretical treatment. To accomplish this we provide a generic result that highlights how well columns of \(A\) selected from an arbitrary matrix with orthonormal rows can be used to approximate \(A\). We then incorporate the results from 3 to provide end to end guarantees for mGKS.

4.1 The mGKSframework↩︎

We briefly summarize mGKSas a framework for selecting columns of a matrix in 1. It is deliberately written to return \(\mathcal{C}\subset [n]\) of cardinality \(k\) representing the selected columns of \(A\). These columns can then be used as needed, e.g., to compute a low-rank approximation.

The rationale for framing mGKSas a framework is that there are lots of choices for how to compute \(\hat{V}_k\). The details of which scheme is used are inconsequential for the theoretical bounds, all that matters is how well \(\hat{V}_k\) can be used to approximate \(A\) (relative to \(V_k\)). Natural choices are to use a randomized SVD [30] (as in [5]) or a sketch of \(A\) whose rows are then orthogonalized (as in [22]); other choices such as Krylov subspace based methods are also good.

4.2 Approximation using an arbitrary subspace↩︎

Given an arbitrary matrix \(Z\in\mathbb{R}^{n\times k}\) (here, one should think of it as an approximation of \(V_k\) as in 1), how well can we approximate \(A\) using columns selected from \(Z^T\)? We address this question in two stages. First, in 2 we prove a generalization of Theorem 9.1 in [30] that allows us to incorporate the overlap between a sampling matrix and an arbitrary orthogonal matrix (as opposed to just the right singular vectors); this result is of independent interest. We then specialize the result in 2 3 to sampling with a subset of the identity matrix (i.e., column selection)—this recovers lemma 4.1 from [31] in the two norm case and yields new bounds in the Frobenius norm.⁶

We now refine 2 when \(\Omega\) is \(k\) columns from the identity matrix, specifically those chosen by a pivoted QR factorization applied to \(Z^T\) denoted \[\label{eqn:ZpQR} Z^T\begin{bmatrix} \Pi_1 & \Pi_2\end{bmatrix} = Q_1 \begin{bmatrix} R_{11} & R_{12}\end{bmatrix}.\tag{8}\] We first consider the two norm and then the Frobenius norm.

[cor:Z2,cor:ZF] show how the overall approximation error can be decomposed into two parts: (1) a part dependent on the approximating subspace and (2) a part due to the column selection scheme applied to \(Z^T\). Methods to build the approximating subspace \(Z\) are well studied (and continue to receive attention), so we do not dwell on them here.

4.3 Theoretical guarantees of mGKS↩︎

In 3 we collect the above results to provide complete control of the approximation efficacy of mGKSin terms of the quality of \(\hat{V}_k\). The theorem deliberately highlights the suboptimality factor of the approximation relative to the SVD. The special case \(\hat{V}_k = V_k\) simplifies the bounds to make them more directly comparable to certain prior work.

While the aforementioned results are quite strong, it is worth highlighting that even when \(\hat{V}_k^T\neq V_k\) is used to select columns \(\mathcal{C}\) the bound \[\|A-P_\mathcal{C}A\|_2 \leq \frac{\sigma_{k+1}(A)}{\sigma_{\min}(V_k(\mathcal{C},:))}\] always holds. While \(\sigma_{\min}(V_k(\mathcal{C},:))\) cannot be easily computed without access to \(V_k\) it is possible this upper bound is smaller than the one from 2. This suggests that 3 could be quite loose in certain situations. The key question to address is when does picking columns from \(\hat{V}_k^T\) ensure good columns of \(V_k^T\) are selected.⁷ Addressing this question directly is quite challenging.

An important aspect of 3 is that it is near optimal in some sense. Specifically, in a slightly limited setting (\(m=n+1\)) [32] shows that there is a matrix \(A\) such that any selection of \(k\) columns yields an approximation at least a factor of \(\sqrt{n/k}\) worse than optimal. Unlike in the Frobenius norm, it is impossible to have a bound independent of \(n\).

4.4 Basis quality↩︎

We have thus far focused on approximation error. When true right singular vectors are used (i.e., \(\hat{V}_k = V_k\)) it is immediate (see, e.g., [1], [4]) that \[\sigma_{\min}(A({:},\mathcal{C}))\geq \sigma_{k}(A)\sigma_{\min}{(V_k(\mathcal{C},{:}))}.\] Using Bischof-Stewart pivoting then implies that \[\sigma_{\min}(A({:},\mathcal{C}))\geq \sigma_{k}(A)/\sqrt{1+k(n-k)}.\]

It is natural to consider what happens when an approximation to \(V_k\) is used. A brief computation shows that for an arbitrary \(\hat{V}_k\) used within mGKSwe have that \[\sigma_{\min}(A({:},\mathcal{C}))\geq \sigma_{\min}{(A\hat{V}_k)}\sigma_{\min}{(\hat{V}_k(\mathcal{C},{:}))} - \eta_2 \sigma_{k+1}(A),\] where \(\eta_2\) is the same as in ?? . Specifically, using Weyl’s inequality we have that \[\begin{align} \sigma_{\min}(A({:},\mathcal{C})) &= \sigma_{\min}(A\hat{V}_k(\hat{V}_k(\mathcal{C},{:}))^T + A(I-\hat{V}_k\hat{V}_k^T)I({:},\mathcal{C})) \\ &\geq \sigma_{\min}(A\hat{V}_k)\sigma_{\min}(\hat{V}_k(\mathcal{C},{:})) - \|A(I-\hat{V}_k\hat{V}_k^T)\|_2 \\ &\geq \sigma_{\min}(A\hat{V}_k)\sigma_{\min}(\hat{V}_k(\mathcal{C},{:})) - \eta_2 \sigma_{k+1}(A). \end{align}\] Notably, as \(\hat{V}_k\) fails to capture the dominant right singular subspace of \(A\) the quality of the basis may degrade accordingly. In extreme cases this bound could even become vacuous, particularly if there is not a gap between \(\sigma_k(A)\) and \(\sigma_{k+1}(A)\).

5 Randomized Bischof-Stewart pivoting for matrices with orthonormal rows↩︎

A byproduct of the proof of 1 is a clear path to a highly efficient randomized variant of Bischof-Stewart pivoting that provides a set of columns with the same theoretical guarantees as the deterministic scheme when applied to matrices with orthonormal rows—not just in expectation, but every time. This method is best suited to the case where \(m\ll n,\) the typical mGKSregime.

5.1 The algorithm↩︎

The key observation from 1 is that to satisfy the overall bounds on growth of \(\|R_{11}^{-1}\|_{\{2,F\}}\) (and the individual singular value bounds) we only need to ensure that at each step \[\label{eqn:RandAccept} f_{i+1} \leq \frac{m-i+1}{m-i}f_i + \frac{n-2i}{m-i},\tag{9}\] where, as before, \(f_i = \|R_{11}^{-1}\|_F^2\) when \(R_{11}\) is \(i\times i\) and \(f_0=0\). When \(A\) has orthonormal rows we are guaranteed an index \(j\) in 5 that satisfies this criterion. However, there are likely many such indices and we need not worry about finding the minimizer. Any index that satisfies the growth bound suffices.

The natural way to leverage this observation is to avoid searching over all columns for a suitable index \(j\) to use as the next pivot and, crucially, to not maintain “state” for all the columns (the expensive part). In the context of 3 this means that we only compute \(R_{11}\) and a representation of \(\Pi\). At any given step only certain columns of \(R_{12}\) and \(M\) are constructed. It is important to note that this means the randomized algorithm, by default, only returns \(\Pi,\) \(Q,\) and \(R_{11}\). If desired, \(Q^T\) can be used to compute \(R_{12}\).

Building off this observation, the randomized algorithm proceeds in the following manner: we randomly sample \(k_b\) columns from \(A,\) update them to the current state using saved Householder reflectors from earlier steps, start selecting columns from the block progressively using Bischof-Stewart pivoting until we are unable to find one that satisfies 9 , sample a new set of \(k_b\) columns and continue till we have selected \(m\) total columns. This process is encoded in 2 and we entitle the algorithm randBSQR. By construction, the selected subset of columns satisfies the same bounds as the deterministic variant; this is encapsulated in 4. For simplicity the algorithm is written to always run for \(m\) steps.

5.2 Experimental setup↩︎

We provide an optimized implementation of 2 as C++ MEX file usable from MATLAB [33].⁸ The implementations, tests, and code to reproduce figures are available at https://github.com/asdamle/BischofStewartQR. All experiments were performed on a MacBook Air with an M3 processor and 16 GB of RAM using the Apple Accelerate framework. Timings are computed using timeit to ensure stable results.

5.3 Matrices↩︎

Throughout these experiments we will primarily use three different test matrices with orthogonal rows. All originate from an \(m\times n\) matrix \(G\) with iid \(\mathcal{N}(0,1)\) entries. The gaussian matrices are constructed by orthogonalizing the rows. To build spiked leverage matrices and needle (simulating a "needle in a haystack" problem) we first weight the columns before orthogonalizing the rows. In the spiked leverage case we place a weight of 100 on \(\lceil 3m/2\rceil\) columns and 1 on the remainder. In the needle example we place a weight of 1 on \(\lceil 5m/4\rceil\) columns and \(10^{-6}\) on all others.

A few more test matrices appear sporadically. The graded leverage matrices start with \(G\) as above and then weight the columns with logarithmically spaced weights between 1 and \(10^{-3}\) before orthogonalizing the rows. The coherent matrices are built by picking \(\lceil m/4 \rceil\) random vectors in \(\mathbb{R}^m\) and then building a matrix with \(n\) columns by uniformly at random sampling one of those vectors to produce each column of \(A,\) and adding \(5\%\) noise to each column and then orthogonalizing the rows. The collinear cluster matrices again start from \(G,\) but replace the first 20% of the columns with a large fixed vector \(v_0\) (chosen at random) plus a small amount of noise, and then orthogonalize the rows. Finally, the chebyshev matrices sample the first \(m\) Chebyshev polynomials at \(n\) random points on \([-1,1]\) and then orthogonalizing the rows.

5.4 Parameter choices↩︎

There are a few parameter choices worth discussing. Selecting \(k_b\) columns at a time yields practical efficiency, as does allowing for the selection of multiple columns out of each block. However, there are some points of tension highlighted by 3. First, we see that there is a trade-off between efficiency and \(\|R_{11}^{-1}\|_F\). While the bound on \(\|R_{11}^{-1}\|_F\) is always satisfied, in some cases there is a small benefit to picking larger blocks. From an efficiency perspective the optimal block size is around \(k_b = 32,\) though that varies slightly with \(m\). If \(k_b\) is too large we end up sampling many more columns in total. Practically, we find \(k_b = m\) to provide a good balance between performance and conditioning. There is also a choice in terms of the sampling distribution for columns. Sampling proportional to the initial squared column norms (i.e., \(\|A({:},j)\|_2^2\)) is a good choice and only adds a one time pre-processing step of \(\mathcal{O}(mn)\). We refer to this scheme as norm-weighted sampling. In certain circumstances uniform sampling could be viable, but in adversarial settings it leads to sampling many more columns than necessary as illustrated in 4.

5.5 Computational cost↩︎

The deterministic variant of Bischof-Stewart pivoting has the same asymptotic cost⁹ as a QR factorization: \(\mathcal{O}(mn\min{(m,n)})\). One drawback of the randomized setting is that we cannot fully leverage the efficient updates of \(R_{11}^{-1}R_{12}\) or column norms of \(M\) since we are not tracking those quantities. Instead, each time we sample a block of columns we must compute its portion of \(R_{12},\) \(R_{11}^{-1}R_{12},\) and \(M\). Assume we have selected \(\ell\) columns so far, then this requires applying \(\ell < m\) Householder reflectors, \(\mathcal{O}(m \ell k_b)\), and running backwards substitution with \(k_b\) right hand sides, \(\mathcal{O}(\ell^2 k_b)\). Running Bischof-Stewart pivoting forward costs at most \(\mathcal{O}(m k_b \min{(m,k_b)})\) (assuming we get \(\min{(m,k_b)}\) pivots out of the sample.) The key point is that each time we sample a block of columns the work done on that block is at most \(\mathcal{O}(m^2 k_b)\).

The overall complexity of 2 is, therefore, entirely dependent on how many blocks we have to sample, \(T_s\), to find \(m\) total columns. Unfortunately, this is not so easy to analyze. For uniform sampling it could be quite bad. However, for norm-weighted sampling we observe that across a broad range of test matrices \(T_s\) is a small constant when \(k_b = m\) (see 4). The cost of randomized Bischof-Stewart pivoting can be cleanly split into two parts: a deterministic part and a non-deterministic part. For norm-weighted sampling the fixed setup cost is \(\mathcal{O}(mn),\) where we compute the squared column norms and build a Fenwick tree [34] to allow for efficient sampling without replacement. The non-deterministic cost¹⁰ is \(\mathcal{O}(T_s m^2 k_b + T_s k_b \log{n}),\) where the \(T_s k_b \log{n}\) arises from sampling the blocks with the aforementioned tree. If \(R_{12}\) is desired, we pick up a \(\mathcal{O}(m^2 n)\) post processing cost.

5.6 Performance↩︎

There are two natural comparison points for randBSQR: (1) how much more efficient is it than the deterministic variant and the standard Golub-Businger pivoting scheme (accessed via DGEQP3), and (2) how does it compare to other randomized methods applicable to this setting such as randomly-pivoted QR (RPQR) [13], [14] (using the efficient implementation from [14] that leverages rejection sampling as in [36], [37]). Notably, RPQR only provides guarantees in expectation on column (or row) based approximation error of an underlying matrix using a GKS like scheme. In fact, in \(\|\cdot\|_F\) it provides stronger guarantees than mGKSowing to its connections to volume sampling. We do not observe such behavior in practice. In contrast, 2 provides guarantees on the quality of the submatrix selected and 4 turns these into approximation error bounds.

Other relevant work includes [38] and [39], neither of which are limited to the case of orthonormal rows. The former guarantees the same output as Golub-Businger QR but avoids excessive updates of columns that will, provably, not be selected. The latter leverages sparse sketches from the right (in contrast to sketches from the left, which are irrelevant when \(m\ll n\) and we want \(m\) columns) to identify a small number of candidate columns from which a good set can then be selected. While similar ideas could, perhaps, be paired with the ideas of 2, the results in [39] suggest the cost of sketching limits the speedup compared to what we observe here. Similarly, while the algorithm in [38] is more general than randBSQR the speedup over DGEQP3 lags what we observe here in the specialized setting.

First, we consider the speedup achieved by 2 relative to both our own optimized deterministic version (see 7.2) and the built-in QR with column pivoting (i.e., DGEQP3 as accessed via MATLAB). 5 shows that for large \(n\) 2 consistently offers a two order of magnitude speedup over the deterministic algorithms when only the selected columns are desired. Even when \(R_{12}\) is desired we still achieve an order of magnitude speedup. This does come at some cost. 6 shows that the quality of the selected subset does lag that of the version that looks at all the columns in some settings. This is expected given the relatively small number of columns that we consider. Nevertheless, the randomized variant always satisfies the theoretical assurances provided.

Now, we explore the relative performance of 2 and RPQR. First, 7 compares the quality of the selected subset. The only theoretical guarantees for RPQR are within the mGKSframework and pertain to the \(\|\cdot\|_F\) approximation accuracy—though we expect a reasonable subset to be found. We observe that randomized Bischof-Stewart consistently outperforms RPQR on metrics related to \(R_{11}^{-1}\). In fact, RPQR regularly fails to satisfy the theoretical bounds 4, which 2 guarantees. Notably, this improved performance does not come at the expense of computational cost. 8 shows that 2 is consistently about twice as fast as RPQR in the asymptotic regime. Because of its connections to volume sampling RPQR benefits from more precise theory about its runtime (in terms of the number of columns that need to be sampled); 2 lacks such assurances. Nevertheless, 8 suggests that \(T_s\) behaves similarly to the sampling requirements for (blocked) RPQR (see [14]).

Finally, we consider how the mGKSframework performs when using randomized Bischof-Stewart or RPQR. To accomplish this we consider synthetic matrices with prescribed spectrum and varying right singular vectors. In 9 we see that even though \(\|R_{11}^{-1}\|\) and, therefore, \(R_{11}^{-1}R_{12}\) behave considerably better for randomized Bischof-Stewart, the projection error \(\|A-P_\mathcal{C}A\|_2\) is often similar for both methods—though they do differ, with randomized Bischof-Stewart outperforming RPQR, around sharp drops in singular values. Randomized Bischof-Stewart does identify a more stable basis for interpolation (as evidenced by the maximum magnitude entry in \(R_{11}^{-1}R_{12}\) being smaller). When an oblique projector is used to determine the low-rank approximation (i.e., up to a permutation from the right our low-rank approximation is \(A({:},\mathcal{C})\begin{bmatrix} I & R_{11}^{-1}R_{12}\end{bmatrix}\)) randBSQR consistently outperforms RPQR.

6 Discussion↩︎

There are several natural questions raised by this work that we briefly discuss. In addition we highlight some applications where the results could have immediate impact.

6.1 Improving Frobenius norm approximation bounds↩︎

Owing to its connections to volume sampling [15]–[17], in the mGKSframework RPQR yields strong theoretical guarantees [13] for Frobenius norm approximation error. A natural question is if Bischof-Stewart satisfies similar (\(n\)-independent) bounds, i.e., if 3 is loose in \(\|\cdot\|_F\). Osinsky shows a greedy strategy can achieve optimal Frobenius norm bounds, but it requires access to the residuals. So, randomness may not be, strictly speaking, necessary. The more pertinent question is if it’s possible to identify a good submatrix of \(V_k\) but observe poor approximation error; while this is not possible in the two norm, the answer is less clear in the Frobenius norm where it is often possible to “hide” error by, e.g., replicating columns.

6.2 Picking rows and columns↩︎

We have framed the most natural application of Bischof-Stewart pivoting as within the mGKSframework, and a clear use of the selected columns is as part of a low-rank approximation. One may also wish to build a two-sided interpolatory low-rank approximation where columns of \(A\) are used to approximate the column space and rows of \(A\) are used to approximate the row space—a \(CUR\) factorization [40]. The path to using Bischof-Stewart in this context is clear. We simply run mGKSon \(A^T\) to identify rows and use those as part of the factorization. Such a scheme readily admits theoretical bounds [31], [40] based on the quality of just a column (or row) based approximation.

6.3 Cholesky factorizations↩︎

There are close connections between pivoted QR factorizations and pivoted Cholesky factorizations. In particular between RPQR and randomly pivoted Cholesky [41], [42]. Therefore, a natural question is what the analogue of Bischof-Stewart pivoting is for Cholesky factorizations—presumably it corresponds to not just pivoting based on (progressively updated) diagonal entries but rather on a combination of diagonal entries and interpolation coefficients.

6.4 Matrices without orthonormal rows↩︎

Applying Bischof-Stewart pivoting to general matrices \(A\) could lead to improved performance over Golub-Businger pivoted QR in certain settings. However, as noted by Stewart [9] the worst case performance is similar. Moreover, the randomized variant cannot be easily extended to this setting as it is unclear what acceptance criterion to use when selecting a column.

There is one setting where Bischof-Stewart is still useful—identifying a well-conditioned square submatrix of a short wide matrix \(W\in\mathbb{R}^{m\times n}\) with \(m < n\). In this case the best possible conditioning of the submatrix is dictated by the singular values of \(W\). While we can no longer prove that Bischof-Stewart can be used to compute a strong rank-revealing factorization, there is a way to identify \(m\) columns of \(W\) whose minimal singular value is as large as we would get from a strong rank-revealing factorization. Specifically, if we compute the reduced QR factorization \(W^T = QR\) and then apply (randomized) Bischof-Stewart pivoting to \(Q^T\) it identifies a subset \(\mathcal{C}\) of columns of \(W\) such that \(\sigma_{\min}(W({:},\mathcal{C})) \geq \sigma_{\min}(W)/\sqrt{1+m(n-m)}\). This result follows immediately from 1 and a singular value variant of Ostrowski’s theorem for eigenvalues [43].

6.5 Applications↩︎

Even with the restriction to orthonormal rows, there are several applications that could immediately benefit from this work. For example, problems in electronic structure theory [44], [45], model reduction [46], and optimal design [47] all directly rely on pivoted QR factorizations of short wide matrices with orthonormal rows. In such applications it is common to appeal to strong rank-revealing QR factorizations for theoretical guarantees, but to simply use DGEQP3 in practice. randBSQR can provide the necessary theoretical assurances while simultaneously providing a means to compute the desired submatrix/columns much more efficiently than is possible with DGEQP3.

7 Bischof and Stewart’s modified pivoting criterion, in detail↩︎

We present the pivoting strategy of Bischof and Stewart in detail. In addition, we provide experimental results from two implementations to show the performance relative to existing, optimized codes (specifically, DGEQP3 leveraging the Apple Accelerate framework). For clarity, we present the algorithm in terms of the linear algebra involved with a focus on how the next pivot is selected—“standard” approaches would be used to ensure computational efficiency. We make a few remarks about the numerical implementation when relevant.

7.1 The algorithm↩︎

Given \(A\in\mathbb{R}^{m\times n}\) for any \(0 \leq i <\min{(m,n)}\) assume we have the partially computed (full) QR factorization \[\label{eq:ApivotedQR} (Q^{(i)})^T A\Pi^{(i)} = \begin{bmatrix}R_{11}^{(i)} & R_{12}^{(i)} \\ & M^{(i)}\end{bmatrix}\tag{10}\] where \(Q^{(i)}\) is orthogonal, \(R_{11}\in\mathbb{R}^{i\times i}\) is upper triangular, \(M\in\mathbb{R}^{(m-i) \times (n-i)}\), and \(R_{12}\in\mathbb{R}^{i\times (n-i)}\). To make this well specified, we adopt the convention that given \(\Pi^{(i)}\) all other matrices are those computed via the first \(i\) steps of a QR factorization without pivoting. By convention \(\Pi^{(0)} = I,\) though the final output of the algorithm is essentially independent of the choice of \(\Pi^{(0)}\) if no ties arise in the pivot selection.

Given \(\Pi^{(i)}\) we describe how to compute \(\Pi^{(i+1)}\). The algorithm then reduces column \(i+1\) of \(A\Pi^{(i+1)}\) to upper triangular form using a Householder reflector, increments \(i\) to \(i+1\) and continues until \(i = \min{(m,n)}\) (unless a partial factorization is explicitly desired). The design choice of Bischof-Stewart is for the first \(i\) columns of \(\Pi^{(i+1)}\) to be those of \(\Pi^{(i)}\) and to choose the \(i+1\) column of \(\Pi^{(i+1)}\) to minimize \(\|(R_{11}^{(i+1)})^{-1}\|_F\) subject to this constraint. 10 provides an algorithmic description of the pivoting scheme; it is written for clarity rather than performance (though it does include the necessary steps to ensure the asymptotic cost is generically \(\mathcal{O}(mn\min{(m,n)})\) and \(\mathcal{O}(mnk)\) if stopped after \(k\) steps). When \(m>n\) the algorithm is easily adapted to return a reduced QR factorization if \(k = \mathop{\mathrm{rank}}{A}\) by dropping rows of \(R\) that are zero.

Diving into 10 we make two key observations. First, as illustrated in 6 the pivoting strategy in [line:BSpiv] due to Stewart [9] explicitly minimizes the growth in the Frobenius norm of \((R_{11}^{(i)})^{-1}\) from one step to the next. Since our convention at the first step is that \(\|w_j\|_2 = 0\), the first pivot selected by this scheme is the same as the first one selected by the classical Golub and Businger scheme. In the unlikely event of a tie at any step any reasonable choice is fine, e.g., picking the smaller index or picking one of the minimizing indices at random.

The second key observation is that [line:wup] is actually maintaining the columns of \((R_{11}^{(i)})^{-1}R_{12}^{(i)}\) (see [9]). This is important to ensure the desired computational scaling—a naı̈ve scheme that recomputes these vectors at each step would pick up an extra factor of \(k\) in the scaling. In fact, the update in [line:wup] can be manipulated to show that we can also cheaply update the squared norm of \(w_{j}\) as \[\label{eq:wnormdowndate} \|\hat{w}_{j}\|_2^2 = \|w_{j}\|_2^2 - 2\alpha_jw_{j}^Tw_{i} + \alpha_j^2(\|w_i\|_2^2+1).\tag{11}\] While useful, this formula is not needed to ensure the correct asymptotic scaling.

7.2 Illustrative performance↩︎

Though it has the same asymptotic complexity, the Bischof-Stewart strategy is more computationally intensive than the method of Golub and Businger (as practically manifest in the LAPACK routine DGEQP3). We provide two optimized implementations of Bischof-Stewart pivoted QR—one C++ mex file usable from MATLAB [33] and one in Julia [49], [50]. The experimental setup mirrors that in 5.

Informally, in the roughly square (\(m\approx n\)) or short and wide setting (\(m \ll n\)) performance can easily be within a factor 2 of the optimized LAPACK routine DGEQP3.¹¹ In fact, in the short wide setting where it is often desirable to compute the interpolation coefficients \(R_{11}^{-1}R_{12}\) the gap is even smaller—a standard pivoted QR requires a subsequent triangular solve to get these coefficients whereas Bischof-Stewart pivoting can get them at minimal additional cost since that quantity is related to the pivoting selection scheme. 11 quantifies these observations in Julia. MATLAB experiments may be found in the main text.

Acknowledgements↩︎

We would like to thank participants of the Simons workshop on complexity and linear algebra for many insightful discussions related to this work. In particular, discussions with Daniel Kressner, Ethan Epperly, Robert Webber, Ilse Ipsen, Mark Embree, Laura Grigori, Yuji Nakatsukasa, Mark Fornace, Alex Townsend, and Michael Lindsey sharpened our thinking of the relationship between mGKSand Osinsky’s results. Implementations were produced using Anthropic’s Claude Code [51]. As was code to run numerical experiments and generate plots. OpenAI Codex [52] was used for code review and early prototyping. Deterministic implementations were validated against a carefully audited oracle implementation (albeit one that was far less performant). Both versions were extensively tested to assess correctness; these tests are included in the repository. AI tools (Anthropic’s Claude and OpenAI’s ChatGPT) were used to aid in proofreading parts of the manuscript and to provide detailed reviews of drafts that guided revisions.

References↩︎