REVIEW 2 major objections 4 minor 52 references
Applied to matrices with orthonormal rows, the Bischof-Stewart pivoting rule provably computes a strong rank-revealing QR factorization, and paired with GKS column selection it matches the guarantees of far more expensive algorithms.
Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →
T0 review · deepseek-v4-flash
2026-08-02 04:56 UTC pith:5A3OT3CW
load-bearing objection Deterministic results are solid — the per-i Bischof-Stewart bounds and mGKS guarantees check out — but the printed randomized algorithm has a broken interpolation-vector update that collapses it to Golub-Businger; fix the pseudocode and this is very publishable. the 2 major comments →
Computing Strong Rank-Revealing Factorizations for Matrices with Orthonormal Rows
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
Core claim
Bischof-Stewart pivoting applied to a matrix with orthonormal rows computes a strong rank-revealing QR factorization: for an n×m Z with Z^T Z = I, the selected columns C satisfy σ_i(Z(C,:)) ≥ 1/√(1 + i(n−m)/(m−i+1)). The proof tracks ||R_11^{-1}||_F^2 through the pivots; because the trailing block M keeps orthonormal rows, the greedy rule becomes a linear recurrence. Fed into GKS column selection, this gives ||A − P_C A||_2 ≤ √(1 + k(n−k)) σ_{k+1}(A) and σ_min(A(:,C)) ≥ σ_k(A)/√(1 + k(n−k)); approximate singular vectors cost a relative factor. A randomized variant, randBSQR, enforces the same per-step bound on sampled blocks and runs about two orders of magnitude faster for m ≪ n.
What carries the argument
The central object is the Bischof-Stewart pivot rule: at each step, choose the column j minimizing (1 + ||R_11^{-1} R_12 e_j||_2^2)/||M e_j||_2^2 — the exact one-step increment, via identity (2.4), of ||R_11^{-1}||_F^2. The load-bearing invariant is that after i pivots of an orthonormal-row Z^T, the trailing block M is (m−i)×(n−i) with orthonormal rows, so ||M||_F^2 = m−i and ||R_11^{-1} R_12||_F^2 = ||R_11^{-1}||_F^2 − i. These two facts turn the greedy minimization into the recurrence f_{i+1} ≤ ((m−i+1)/(m−i)) f_i + (n−2i)/(m−i), whose solution yields the singular-value bounds; the randomized algorithm keeps only this recurrence as an acceptance test, sampling blocks of columns instead of
Load-bearing premise
The load-bearing premise, used in the proof of Theorem 3.1, is the exact-arithmetic invariant that after every pivot the unreduced block M retains orthonormal rows, giving ||M||_F^2 = m−i and ||R_11^{-1} R_12||_F^2 = ||R_11^{-1}||_F^2 − i; floating-point Householder updates can break it, and the randomized acceptance test (5.1) is shielded only by the heuristic slack factor of Remark 4.
What would settle it
Run the deterministic or randomized Bischof-Stewart algorithm in double precision on an orthonormal-row matrix engineered to saturate the bound (e.g., n = m + 1 with a highly correlated, ill-conditioned column structure), record f_i = ||R_11^{-1}||_F^2 after each pivot, and check whether rounding ever pushes the computed σ_min of the selected submatrix below 1/√(1 + i(n−m)/(m−i+1)); a single violation after the factorization completes would show the exact-arithmetic proof does not survive finite precision. Separately, constructing a family of matrices where norm-weighted sampling forces many b
If this is right
- GKS column selection paired with Bischof-Stewart pivoting achieves the same two-norm approximation and basis-conditioning guarantees as applying a strong rank-revealing QR factorization directly to A, at ordinary pivoted-QR cost.
- The per-step singular-value bounds are, for i < m, stronger than those offered by the standard polynomial-time strong rank-revealing algorithm, and no column-swap refinement is required.
- The mGKS guarantees are deterministic — they hold every time — unlike volume-sampling and randomly pivoted schemes whose bounds hold only in expectation.
- When only an approximation V̂_k of the leading right singular vectors is available, the bounds degrade by exactly the relative factor ||A(I − V̂_k V̂_k^T)||/||A(I − V_k V_k^T)||, which makes the framework usable when A is accessible only through matrix-vector products.
- The randomized variant randBSQR returns a subset with the same per-step guarantees, runs up to two orders of magnitude faster than the deterministic scheme for m ≪ n, and is about twice as fast as randomly-pivoted QR in the asymptotic regime.
Where Pith is reading between the lines
- The proof actually licenses a family of 'lazy' pivoting schemes: since any pivot satisfying the per-step recurrence (5.1) inherits the full guarantee, one could accept columns from small candidate pools, partial norm updates, or other cheap surrogates without losing the bounds — an extension the paper gestures at but does not develop.
- The deterministic (not merely in-expectation) character of randBSQR is arguably the sharpest break with volume-sampling methods; replacing the heuristic slack factor of Remark 4 with a backward error analysis would carry the guarantees into finite precision. This is an editorial extrapolation, not a paper claim.
- The paper leaves the Frobenius-norm gap open: volume-sampling schemes can beat the mGKS Frobenius bound, so a natural test is whether Bischof-Stewart's interpolation coefficients R_11^{-1} R_12 carry hidden volume-like structure, as the maximum-volume connection hints.
- A Cholesky analogue — pivoting on a blend of diagonal entries and interpolation coefficients — would supply a deterministic counterpart to randomly pivoted Cholesky for kernel matrices; the paper raises the question without answering it.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper studies the Bischof–Stewart pivoting strategy for matrices with orthonormal rows and shows that it computes a strong rank-revealing QR factorization (Theorem 3.1). It then uses this result in the mGKS column-selection framework to obtain two-norm and Frobenius-norm bounds for column subset selection, including the case where only an approximation of the right singular vectors is available (Theorem 4.4 and Section 4.4). A randomized variant, randBSQR, is proposed (Algorithm 5.1), with the claim that it satisfies the same theoretical guarantees while being much faster; numerical experiments compare it with deterministic Bischof–Stewart pivoting, LAPACK's DGEQP3, and RPQR.
Significance. If the algorithmic claims are correct, the paper makes a strong contribution: it gives a parameter-free, polynomial-cost pivoting scheme with per-step singular-value bounds for orthonormal-row matrices, and it extends GKS to match the two-norm guarantees of strong RRQR at essentially pivoted-QR cost. The randomized variant, if implemented as intended, would convert an expensive global search into a block-sampled procedure with deterministic (not expectation-only) guarantees. The proof of Theorem 3.1 is a clean induction, and the mGKS corollaries follow in a straightforward way from standard CS-decomposition identities. The experimental section is reproducible and the comparison with RPQR is relevant. However, the algorithms as printed do not implement the pivot rule that the theorems analyze, so the paper's central computational claims are not supported without correction.
major comments (2)
- [Appendix A, Algorithm A.1; Section 5.1, Algorithm 5.1] The printed update for w_j is not the Bischof–Stewart pivot rule used in Theorem 3.1. The quantity required in (2.3) is w_j = R11^{-1} R12 e_j. Algorithm A.1 initializes w_j = [] and, after the first Householder step, sets w_j = (w_j - alpha_j w_i)/alpha_j (line 9, alpha_j = R(i,j)/R(i,i)). Since w_j and w_i are empty at i=1, w_j remains empty; from step 2 onward the selection criterion reduces to maximizing ||M e_j||^2, i.e., the Golub–Businger rule. The same defect appears in Algorithm 5.1: line 9 sets w_j = X(1:i,j), which is R12 e_j rather than R11^{-1} R12 e_j, and line 25 applies the same invalid division. The correct block-inverse update is w_j <- [w_j - alpha_j w_ell; alpha_j], appending the new coordinate; the printed formula neither appends the new coordinate nor preserves vector length, and it can divide by zero when R(i,j)=0. Consequently, the quantity c in Algorithm 5.1, lin
- [Section 5.1, Remark 4 and Theorem 5.1] The theoretical guarantees are exact-arithmetic statements, but the paper repeatedly describes Algorithm 5.1 as 'always' satisfying the bounds in floating point (e.g., Section 5.6). The proof of Theorem 3.1 relies on the invariants ||M||_F^2 = m-i and ||R11^{-1}R12||_F^2 = f_i - i after each step; in finite precision, Householder updates and the norm updates in Algorithm 5.1 can break these invariants. The slack factor in Remark 4 is heuristic, not a perturbation theorem. To make the claim that the implemented randomized algorithm satisfies (5.1) rigorous, the paper should either state explicitly that all guarantees are exact-arithmetic and that the numerical experiments are empirical, or provide a backward-error analysis for the acceptance test.
minor comments (4)
- [Section 3, after Theorem 3.1] The extraction of the entrywise bound on R11^{-1}R12 from the proof of Theorem 3.1 is informal. A one-line derivation from f_i - i to the max-norm bound would make the 'strong rank-revealing' statement easier to verify.
- [Algorithm A.1, line 9] Aside from the incorrect update, alpha_j = R(i,j)/R(i,i) can be zero for a remaining column whose residual is nonzero, so the division in the printed formula is not merely an algebraic typo but can cause a runtime failure. The corrected recurrence avoids this issue.
- [Section 5.5] The complexity statement depends on the empirical observation that T_s is a small constant for norm-weighted sampling. The paper acknowledges the worst-case T_s <= m n / k_b, but the 'two orders of magnitude' claim in the abstract and Section 5.6 should be labeled as empirical performance, not as a theorem.
- [Figure 6] The caption contains a formatting issue ('Figure6' without a space) and the figure reference in Section 5.5 should be checked.
Circularity Check
No significant circularity: the core bounds are proved from the inverse-update recurrence (2.4) and standard CS/interlacing arguments; self-citations are background, not load-bearing.
full rationale
Theorem 3.1 is derived rather than assumed: starting from the explicit formula for R_hat_11^{-1}, the paper obtains (2.4), and the invariants ||M||_F^2 = m-i and ||R11^{-1}R12||_F^2 = ||R11^{-1}||_F^2 - i follow from the orthonormal-row structure at each QR step. The inductive bound on f_i and conversion s_i <= f_i - (i-1) are direct algebra; no fitted parameter or cited theorem supplies the result. Corollaries 4.2 and 4.3 are proved in place via a CS decomposition, and Theorem 4.4 only inserts the Theorem 3.1 bound into them. Theorem 5.1 is a correctness theorem for an algorithm that explicitly enforces condition (5.1); the text itself says "By construction, the selected subset satisfies the same bounds," so this is a design guarantee rather than a prediction extracted from fitted data. The floating-point caveat in Remark 4 and the appendix update issue raised in the skeptic note are numerical-correctness concerns, not circularity. Self-citations [6] and [20] are used for background/context or to name a property already proved directly in Theorem 3.1 (with the ||R11^{-1}R12||_max bound extracted from the proof), so they are not load-bearing. Osinsky overlap is disclosed. Overall: no significant circularity; score reflects only minor non-load-bearing self-citations.
Axiom & Free-Parameter Ledger
free parameters (3)
- block size k_b =
default m (e.g. 256 in experiments)
- sampling weights g_j =
default ||A(:,j)||_2^2 / m
- slack factor in acceptance test =
small multiple of machine precision
axioms (6)
- domain assumption Exact arithmetic preserves orthonormal rows of the unreduced block M throughout the partial QR factorization.
- standard math The CS decomposition of [Z, Z_perp]^T [Pi1, Pi2] gives ||W2 W1^-1|| = ||R11^-1 R12|| in both the 2- and Frobenius norms.
- standard math Weyl's inequality and singular-value interlacing justify the basis-quality bound in Section 4.4.
- domain assumption In mGKS, V_hat_k is exactly orthonormal and k < rank(A); the residual ratios eta_2 and eta_F are finite.
- standard math At each step of randomized Bischof-Stewart, at least one remaining column satisfies the growth bound (5.1).
- standard math Ostrowski's theorem connects sigma_min(W(:,C)) to sigma_min(W) after QR reduction in the Section 6.4 generalization.
read the original abstract
We show that a pivoting strategy due to Stewart (based on work by Bischof) computes a strong rank-revealing factorization when applied to a matrix with orthonormal rows. When paired with the classical column selection algorithm of Golub, Klema, and Stewart (GKS) it helps achieve rank-$k$ approximation accuracy bounds and basis conditioning as good as those from applying a strong rank-revealing factorization directly to A. We then extend this framework in two directions: (1) providing analysis of GKS when only approximations of right singular vectors are available and (2) providing a randomized variant of the pivoting strategy for matrices with orthonormal rows that achieves the same theoretical guarantees but can return the desired subset two orders of magnitude faster than the deterministic variant.
Reference graph
Works this paper leans on
-
[1]
https://claude.com/product/claude-code
Anthropic C laude C ode . https://claude.com/product/claude-code
-
[2]
https://julialang.org
Julia . https://julialang.org
-
[3]
https://www.mathworks.com/
MATLAB . https://www.mathworks.com/
-
[4]
https://openai.com/codex/
Open AI C odex . https://openai.com/codex/
-
[5]
N. Amsel, Y. Baumann, P. Beckman, P. Bürgisser, C. Camaño, T. Chen, E. Chow, A. Damle, M. Derezinski, M. Embree, E. N. Epperly, R. Falgout, M. Fornace, A. Greenbaum, C. Greif, D. Halikias, Z. Huang, E. Jarlebring, Y. Koutis, D. Kressner, R. Kyng, J. Liesen, J. Lok, R. A. Meyer, Y. Nakatsukasa, K. Pearce, R. Peng, D. Persson, E. Rebrova, R. Schneider, R. S...
arXiv 2026
-
[6]
Armstrong, A
R. Armstrong, A. Buzali, and A. Damle , Structure-aware analyses and algorithms for interpolative decompositions , SIAM Journal on Scientific Computing, 47 (2025), pp. A1527--A1554
2025
-
[7]
R. Armstrong and A. Damle , Collect, commit, expand: Efficient CPQR -based column selection for extremely wide matrices , arXiv preprint arXiv:2501.18035, (2025)
Pith/arXiv arXiv 2025
-
[8]
Avron and C
H. Avron and C. Boutsidis , Faster subset selection for matrices and applications , SIAM Journal on Matrix Analysis and Applications, 34 (2013), pp. 1464--1499
2013
-
[9]
Barthelm\'e, N
S. Barthelm\'e, N. Tremblay, and P.-O. Amblard , A faster sampler for discrete determinantal point processes , in Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, F. Ruiz, J. Dy, and J.-W. van de Meent, eds., vol. 206 of Proceedings of M achine L earning R esearch, PMLR, 25--27 Apr 2023, pp. 5582--5592
2023
-
[10]
Belhadji, R
A. Belhadji, R. Bardenet, and P. Chainais , A determinantal point process for column subset selection , Journal of machine learning research, 21 (2020), pp. 1--62
2020
-
[11]
Bezanson, A
J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah , Julia: A fresh approach to numerical computing , SIAM R eview, 59 (2017), pp. 65--98
2017
-
[12]
C. H. Bischof , Incremental condition estimation , SIAM Journal on Matrix Analysis and Applications, 11 (1990), pp. 312--322
1990
-
[13]
Boutsidis, P
C. Boutsidis, P. Drineas, and M. Magdon-Ismail , Near-optimal column-based matrix reconstruction , SIAM Journal on Computing, 43 (2014), pp. 687--717
2014
-
[14]
Businger and G
P. Businger and G. H. Golub , Linear least squares solutions by H ouseholder transformations , Numerische Mathematik, 7 (1965), pp. 269 -- 276
1965
-
[15]
Chandrasekaran and I
S. Chandrasekaran and I. C. F. Ipsen , On rank-revealing factorisations , SIAM Journal on Matrix Analysis and Applications, 15 (1994), pp. 592--622
1994
-
[16]
L. Chen, Z. Liu, W. He, and B. Dong , Iteris: Agentic research loops for computational mathematics , arXiv preprint arXiv:2606.02484, (2026)
Pith/arXiv arXiv 2026
-
[17]
Y. Chen, E. N. Epperly, J. A. Tropp, and R. J. Webber , Randomly pivoted C holesky: Practical approximation of a kernel matrix with few entry evaluations , Communications on Pure and Applied Mathematics, 78 (2025), pp. 995--1041
2025
-
[18]
Cortinovis and D
A. Cortinovis and D. Kressner , Adaptive randomized pivoting for column subset selection, DEIM , and low-rank approximation , SIAM Journal on Matrix Analysis and Applications, 47 (2026), pp. 25--47
2026
-
[19]
Cortinovis and L
A. Cortinovis and L. Ying , A sublinear-time randomized algorithm for column and row subset selection based on strong rank-revealing QR factorizations , SIAM Journal on Matrix Analysis and Applications, 46 (2025), pp. 22--44
2025
-
[20]
Damle, S
A. Damle, S. Glas, A. Townsend, and A. Yu , Estimating a matrix's singular values with interpolative decompositions , Linear Algebra and its Applications, 731 (2026), pp. 306--342
2026
-
[21]
Damle and L
A. Damle and L. Lin , Disentanglement via entanglement: a unified method for W annier localization , Multiscale Modeling & Simulation, 16 (2018), pp. 1392--1410
2018
-
[22]
Damle, L
A. Damle, L. Lin, and L. Ying , Compressed representation of K ohn-- S ham orbitals via selected columns of the density matrix , Journal of chemical theory and computation, 11 (2015), pp. 1463--1469
2015
-
[23]
Derezi \'n ski, K
M. Derezi \'n ski, K. L. Clarkson, M. W. Mahoney, and M. K. Warmuth , Minimax experimental design: Bridging the gap between statistical and worst-case approaches to least squares regression , in Conference on Learning Theory, PMLR, 2019, pp. 1050--1069
2019
-
[24]
Deshpande and L
A. Deshpande and L. Rademacher , Efficient volume sampling for row/column subset selection , in 2010 IEEE 51st Annual Symposium on Foundations of Computer Science, 2010, pp. 329--338
2010
-
[25]
Deshpande, L
A. Deshpande, L. Rademacher, S. S. Vempala, and G. Wang , Matrix approximation and projective clustering via volume sampling , vol. 2, Theory of Computing, 2006, pp. 225--247
2006
-
[26]
Deshpande and S
A. Deshpande and S. Vempala , Adaptive sampling and fast low-rank matrix approximation , in International Workshop on Approximation Algorithms for Combinatorial Optimization, Springer, 2006, pp. 292--303
2006
-
[27]
Y. Dong, C. Chen, P.-G. Martinsson, and K. Pearce , Robust blockwise random pivoting: Fast and accurate adaptive interpolative decomposition , SIAM Journal on Matrix Analysis and Applications, 46 (2025), pp. 1791--1815
2025
-
[28]
Dong and P.-G
Y. Dong and P.-G. Martinsson , Simpler is better: a comparative study of randomized pivoting algorithms for CUR and interpolative decompositions , Advances in Computational Mathematics, 49 (2023), p. 66
2023
-
[29]
Drineas and I
P. Drineas and I. C. Ipsen , Low-rank matrix approximations do not need a singular value gap , SIAM Journal on Matrix Analysis and Applications, 40 (2019), pp. 299--319
2019
-
[30]
Drma c and Z
Z. Drma c and Z. Bujanovi \'c , On the failure of rank-revealing QR factorization software--a case study , ACM Transactions on Mathematical Software (TOMS), 35 (2008), pp. 1--28
2008
-
[31]
Drma c and S
Z. Drma c and S. Gugercin , A new selection operator for the discrete empirical interpolation method---improved a priori error bound and extensions , SIAM Journal on Scientific Computing, 38 (2016), pp. A631--A648
2016
-
[32]
E. N. Epperly , Adaptive randomized pivoting and volume sampling , arXiv preprint arXiv:2510.02513, (2025)
Pith/arXiv arXiv 2025
-
[33]
E. N. Epperly, J. A. Tropp, and R. J. Webber , Embrace rejection: Kernel matrix approximation by accelerated randomly pivoted cholesky , SIAM Journal on Matrix Analysis and Applications, 46 (2025), pp. 2527--2557
2025
-
[34]
Eswar, V
S. Eswar, V. Rao, and A. K. Saibaba , Bayesian d-optimal experimental designs via column subset selection , SIAM Journal on Scientific Computing, 48 (2026), pp. A905--A928
2026
-
[35]
I. Fakih and L. Grigori , Efficient QR -based column subset selection through randomized sparse embeddings , arXiv preprint arXiv:2509.03198, (2025)
Pith/arXiv arXiv 2025
-
[36]
P. M. Fenwick , A new data structure for cumulative frequency tables , Software: Practice and Experience, 24 (1994), pp. 327--336
1994
-
[37]
Golub , Numerical methods for solving linear least squares problems , Numerische mathematik, 7 (1965), pp
G. Golub , Numerical methods for solving linear least squares problems , Numerische mathematik, 7 (1965), pp. 206--216
1965
-
[38]
G. H. Golub, V. Klema, and G. Stewart , Rank degeneracy and least squares problems , Tech. Rep. STAN-CS-76-559, Stanford University, 1976
1976
-
[39]
S. A. Goreinov, E. E. Tyrtyshnikov, and N. L. Zamarashkin , A theory of pseudoskeleton approximations , Linear algebra and its applications, 261 (1997), pp. 1--21
1997
-
[40]
L. Grigori and Z. Xue , Randomized strong rank-revealing QR for column subset selection and low-rank matrix approximation , arXiv preprint arXiv:2503.18496, (2025)
Pith/arXiv arXiv 2025
-
[41]
Gu and S
M. Gu and S. C. Eisenstat , Efficient algorithms for computing a strong rank-revealing QR factorization , SIAM Journal on Scientific Computing, 17 (1996), pp. 848--869
1996
-
[42]
Guruswami and A
V. Guruswami and A. K. Sinop , Optimal column-based low-rank matrix reconstruction , in Proceedings of the twenty-third annual ACM-SIAM symposium on Discrete Algorithms, SIAM, 2012, pp. 1207--1214
2012
-
[43]
Halko, P
N. Halko, P. G. Martinsson, and J. A. Tropp , Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions , SIAM Review, 53 (2011), pp. 217--288
2011
-
[44]
Hong and C.-T
Y. Hong and C.-T. Pan , Rank-revealing QR factorizations and the singular value decomposition , Mathematics of Computation, 58 (1992), pp. 213 -- 232
1992
-
[45]
Kahan , Numerical linear algebra , Canadian Mathematical Bulletin, 9 (1966), pp
W. Kahan , Numerical linear algebra , Canadian Mathematical Bulletin, 9 (1966), pp. 757--801
1966
-
[46]
Kozyrev and A
I. Kozyrev and A. Osinsky , Subset selection for matrices in spectral norm , BIT Numerical Mathematics, 66 (2026), p. 45
2026
-
[47]
Liberty, F
E. Liberty, F. Woolfe, P.-G. Martinsson, V. Rokhlin, and M. Tygert , Randomized algorithms for the low-rank approximation of matrices , Proceedings of the National Academy of Sciences, 104 (2007), pp. 20167--20172
2007
-
[48]
M. W. Mahoney and P. Drineas , CUR matrix decompositions for improved data analysis , Proceedings of the National Academy of Sciences, 106 (2009), pp. 697--702
2009
-
[49]
Osinsky , Close to optimal column approximation using a single SVD , Linear Algebra and its Applications, 725 (2025), pp
A. Osinsky , Close to optimal column approximation using a single SVD , Linear Algebra and its Applications, 725 (2025), pp. 359--377
2025
-
[50]
A. M. Ostrowski , A quantitative formulation of sylvester's law of inertia , Proceedings of the National Academy of Sciences, 45 (1959), pp. 740--744
1959
-
[51]
D. C. Sorensen and M. Embree , A DEIM induced CUR factorization , SIAM Journal on Scientific Computing, 38 (2016), pp. A1454--A1482
2016
-
[52]
G. W. Stewart , Incremental condition calculation and column selection , Tech. Rep. UMIACS TR 90-87, CS TR 2495, University of Maryland, 1990
1990
discussion (0)
Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.