{"id":"ab567785-9b35-4962-b332-955a9c23371d","arxiv_id":"2508.20269","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"low","formal_verification":"none","parameter_count":2,"one_line_summary":"A randomized Golub-Kahan algorithm and new randomized LSQR, CGLS, and LSMR solvers, with hybrid Tikhonov variants, are proposed for rectangular linear inverse problems.","lead":"This paper develops randomized versions of Golub-Kahan and related Krylov solvers that replace expensive inner products with cheap sketched approximations for large-scale inverse problems. The new methods work for rectangular matrices and include automatic regularization-parameter selection.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The numerical accuracy claim depends on the stability of an un-reorthogonalized, sketched classical Gram-Schmidt step that the paper supports only with an undocumented empirical remark; even with large sketches, loss of Θ-orthogonality can invalidate the optimality equations (2.20)-(2.26).","rationale":"The exact-arithmetic derivation of the rGK factorizations and the resulting solvers is internally consistent, and the residual-optimality arguments are correct under the stated orthonormality assumptions. The reader's weakest assumption about sketch dimensions is relevant, but I see a more specific load-bearing condition: the numerical reliability of the simplified, un-reorthogonalized randomized Gram-Schmidt procedure used in Algorithm 2.2. Even if ℓm and ℓn are large enough to make the sketches good subspace embeddings, classical Gram-Schmidt in floating-point arithmetic can lose orthogonality, and the paper's only support for using the simplified version is an undocumented empirical claim in Section 2.1. The paper also admits in Section 5 that a full complexity analysis and highly efficient implementation are still in progress, which strengthens the case for a conditional rather than unconditional acceptance. The reader's concerns about single runs, missing runtime experiments, and code availability are valid but secondary; the stability of the orthogonalization is more directly tied to whether the central factorization-based optimality claims hold for the actual computed iterates. The proposed test would settle this by measuring orthogonality loss and factorization residuals on a realistic ill-conditioned problem. Since the issue is addressable and the exact-arithmetic framework is sound, the appropriate verdict remains CONDITIONAL, so I do not change the reader's verdict.","tokens_in":21067,"tokens_out":17124,"duration_ms":163874,"concrete_test":"On the seismic tomography example (m=131072, n=65536, K=100), run Algorithm 2.2 with a large sketch, e.g., ℓm=5243 and ℓn=2622, and at every iteration compute the loss of Θ-orthogonality E_U(k)=||S_{k+1}^T S_{k+1}-I||_F and E_V(k)=||P_{k+1}^T P_{k+1}-I||_F, where S_{k+1}=Θ(m)U_{k+1} and P_{k+1}=Θ(n)V_{k+1}. Also compute the factorization residuals ||AV_k-U_{k+1}M_k||_F/||A||_F and ||A^T U_{k+1}-V_{k+1}T_{k+1}||_F/||A||_F, and compare them against the same quantities for a variant with one step of reorthogonalization. If the orthogonality loss and factorization residuals stay near machine precision, the stability concern does not land; if they grow with k, the simplified Gram-Schmidt step is the weak point and the solvers need reorthogonalization or a more stable coefficient solve.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The derivations of rLSQR, rCGLS, and rLSMR and their hybrid variants rest on the rGK factorizations AV_k = U_{k+1}M_k and A^T U_{k+1} = V_{k+1}T_{k+1} in (2.14). These factorizations are exact only if the columns of Θ(m)U_{k+1} and Θ(n)V_{k+1} are orthonormal, and Algorithm 2.2 produces them by a single, un-reorthogonalized classical Gram-Schmidt sweep with coefficients computed only from sketched inner products (steps 13-15 and 21-23). The paper explicitly acknowledges in Section 2.1 that more stable implementations exist, but justifies the simpler version with 'we did not observe substantial differences when testing on a variety of inverse problems' without presenting those tests. In floating-point arithmetic, classical Gram-Schmidt loses orthogonality at a rate governed by the conditioning of the basis, and for severely ill-conditioned inverse problems this can make the factorizations (2.14) inaccurate. When that happens, the identities used in (2.20), (2.24), and (2.26) no longer describe the actual iterates, so the residual-optimality arguments collapse. The erratic behavior at small sketch sizes (Figures 5 and 7) is consistent with loss of Θ-orthogonality, but the paper never measures orthogonality loss, so it is impossible to tell whether the failures come from insufficient sketching or from unstable orthogonalization. This is a load-bearing gap in the central accuracy claim.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper develops randomized Krylov subspace methods for large-scale linear inverse problems with rectangular forward operators A. It extends the randomized Gram-Schmidt approach of Balabanov and Grigori to a randomized Golub-Kahan (rGK) algorithm, which uses two sketching matrices and produces the factorizations AV_k = U_{k+1}M_k and A^T U_{k+1} = V_{k+1}T_{k+1}. On top of these factorizations the authors derive randomized LSQR, CGLS, and LSMR solvers and hybrid Tikhonov-regularized variants with automatic parameter selection via discrepancy principle and weighted GCV. The paper also reviews the randomized Arnoldi/GMRES setting and presents numerical experiments for image deblurring and seismic tomography.","tokens_in":21384,"tokens_out":5412,"duration_ms":50793,"significance":"If the algorithmic claims hold, this is a useful contribution: it extends randomized inner-product-free Krylov techniques from square to rectangular inverse problems, provides new projected solvers with exact-arithmetic derivations, and introduces hybrid variants with adaptive regularization. The derivations of the factorizations and projected systems in Section 2 are careful and self-contained, and the paper is honest about its scope, including the lack of a theoretical analysis of the low-rank approximations. The numerical experiments do show that the methods can work when the sketch sizes are sufficiently large. The main weakness is that the central accuracy and optimality claims depend on the preservation of sketched orthonormality by an un-reorthogonalized classical Gram-Schmidt process, and this is not analyzed or tested. The efficiency motivation is also supported only by flop counts, not by runtime experiments.","major_comments":[{"comment":"The derivations of rLSQR, rCGLS, and rLSMR are exact only if the columns of Θ(m)U_{k+1} and Θ(n)V_{k+1} are orthonormal, since identities such as (2.20) and (2.27) use that orthonormality explicitly. Algorithm 2.2 computes these bases by a single classical Gram-Schmidt sweep with coefficients and norms obtained from sketched vectors (steps 13-17 and 21-25). In floating-point arithmetic, classical Gram-Schmidt can lose orthogonality at a rate governed by the conditioning of the basis, and for severely ill-conditioned inverse problems this can make the factorizations (2.14) inaccurate. The only evidence offered for the stability of the simplified orthogonalization is the undocumented empirical remark in Section 2.1 that no substantial differences were observed on a variety of inverse problems. Since the residual-optimality properties of the proposed solvers collapse when the sketched orthonormality is lost, the paper should either add a reorthogonalized or otherwise stabilized variant, report quantitative measurements of the loss of Θ-orthogonality, or provide an error analysis that bounds the effect of this loss on the factorizations and projected solutions. The erratic behavior at small sketch sizes in Figures 5 and 7 is consistent with this concern, but the paper does not distinguish it from insufficient sketching.","section":"Section 2.2, Algorithm 2.2, Eqs. (2.14), (2.19)-(2.27)"},{"comment":"The claimed computational advantage of rGK over GKB and reorthogonalized GKB is based entirely on flop-count estimates and the synthetic cost comparisons in Figure 1. The numerical experiments in Section 4 report only relative reconstruction errors and do not include any wall-clock times, iteration timings, or comparisons of actual runtime for rGK versus GKB/ro-GKB. Because the main motivation stated in the introduction is the reduction of computational cost, the paper needs at least a basic runtime validation on the reported test problems to confirm that the flop-count advantages translate into practice, especially given the full reorthogonalization and storage requirements of the randomized bases.","section":"Section 4, cost claims in Section 2.2, Eq. (2.16)"},{"comment":"The discrepancy principle and weighted GCV are adapted to randomized hybrid methods by replacing the true residual with the sketched residual Θ(m)(b - Ax_k(λ)). No error analysis is provided that relates the λ selected by these sketched criteria to the λ that would be selected with the true residual, and the numerical results in Figure 11 show that the GCV-based parameter can visibly overestimate the optimal value. Since the adaptive selection of λ is one of the advertised features of the hybrid methods, the heuristic nature of these modifications should be stated more prominently, or a perturbation bound in terms of the embedding quality ε and the sketch dimensions should be supplied.","section":"Section 3, Eqs. (3.15)-(3.16)"}],"minor_comments":[{"comment":"The numerator of the WGCV functional is written as n ||Θ(m)(b - Ax_k(λ))||^2, but the standard GCV formula for an m-dimensional data space uses the factor m; since the denominator traces involve the m×m matrix AA_reg^†(λ), the factor n appears to be a typo and should be m or justified explicitly.","section":"Section 3, Eq. (3.16)"},{"comment":"The sentence 'Similar results are observed for rGMRES-GCV (shown in the right plot of Figure 4)' should refer to Figure 5, since Figure 4 is the comparison with deterministic GMRES and Figure 5 is the one with varying ℓ_n for hybrid rGMRES with GCV.","section":"Section 4.1, Figure 5 and surrounding text"},{"comment":"There is a typo in 'we use the SHRT' which should read 'SRHT'.","section":"Section 4.1"},{"comment":"The notation e_{k+1}^T for the k×(k+1) matrix [I_k, 0]T_{k+1} is confusing, because e_{k+1} is elsewhere used for the (k+1)-st canonical basis vector; a different symbol, for instance T_{k+1}^{(k)}, would improve readability.","section":"Section 2.3, Eqs. (2.24)-(2.26)"}],"recommendation":"major_revision","confidential_remarks":"The paper is a solid algorithmic extension of randomized Gram-Schmidt ideas to rectangular inverse problems, and the exact-arithmetic derivations appear correct. The main risk is the finite-precision stability of the proposed rGK factorization: the paper's central optimality claims rely on sketched orthonormality that is produced by an un-reorthogonalized classical Gram-Schmidt pass, and the supporting numerical evidence is only anecdotal. In my view this is fixable within the manuscript's scope by adding orthogonality measurements, a reorthogonalized variant, or a stability analysis, together with runtime experiments and a small set of modifications to the GCV discussion. I therefore recommend major revision rather than rejection."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Quick take: this is a real contribution, not a repackaging. The randomized Golub-Kahan factorization for rectangular A and the rLSQR, rCGLS, and rLSMR solvers built on it are new algorithmic constructions; the hybrid Tikhonov variants are a natural but useful extension. The factorizations (2.14) are derived carefully and the optimality arguments for rLSQR and rLSMR check out in exact arithmetic, assuming the sketched bases are Θ-orthonormal. The experiments on 262k-pixel deblurring and 131k x 65k seismic tomography support the claimed behavior, and it is to the authors' credit that they include a 24-run box plot for one experiment and openly say in Section 5 that the efficient implementation and complexity analysis are still in progress.\n\nWhere it gets soft: the stress-test note lands. Algorithm 2.2 uses a single un-reorthogonalized classical Gram-Schmidt sweep with sketched inner products, and the paper justifies this with the remark in Section 2.1 that more stable implementations did not produce substantial differences 'when testing on a variety of inverse problems' - without showing those tests. That is a missing piece of evidence, not a fatal flaw. The erratic small-sketch behavior in Figures 5 and 7 is consistent with loss of Θ-orthogonality, but the paper never measures the orthogonality loss, so the reader cannot tell whether the failure is from an undersized sketch or from unstable orthogonalization. The central accuracy claim is conditional on that. A second soft spot is the flop-count story: the cost estimates are plausible, but there are no runtime experiments and no communication-protocol measurement, so the headline 'inner-product bottleneck' benefit is asserted, not demonstrated. Minor: the sketched GCV and discrepancy principle are applied without error analysis or convergence guarantees, and the code is promised but not yet available. The random runs are mostly single trials; the one box-plot experiment shows spread, so the reader should treat the error curves as indicative.\n\nNet: the paper does what it claims - it extends randomized Gram-Schmidt to rectangular inverse problems and delivers solvers that behave like their deterministic counterparts when the sketch is large enough. The stability gap means the theory is not as tight as the presentation implies, but it is addressable with reorthogonalization or at least diagnostic measurements.\n\nWho benefits: computational inverse problems researchers and anyone using Krylov methods where inner products are costly. I would cite it and would bring it to reading group. Recommendation: send to referees. Ask for orthogonality diagnostics, runtime comparisons, code release, and ideally a second look at the sketched parameter-choice rules. This is a conditional accept, not a reject.","headline":"A solid methods paper: the rGK factorization and its three solvers are genuinely new, the derivations hold up in exact arithmetic, and the main gaps are unquantified stability of the sketched Gram-Schmidt and missing runtime evidence.","tokens_in":21926,"tokens_out":2996,"would_cite":true,"duration_ms":26769,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65F22","65F10","65K10","15A29"],"pacs":[],"model":"deepseek-v4-flash","headline":"Randomized Golub-Kahan bidiagonalization extends randomized Krylov solvers to rectangular inverse problems, producing regularized solutions comparable to deterministic Krylov methods when the sketch dimensions are large enough.","keywords":["randomized Golub-Kahan","randomized Krylov methods","hybrid projection methods","Tikhonov regularization","inverse problems","randomized Gram-Schmidt","LSQR","LSMR"],"falsifier":"Run Algorithm 2.2 on a rectangular inverse problem, such as the seismic tomography example with $A$ of size $131072\\times 65536$, while increasing the sketch dimensions $\\ell_m/\\ell_n$ across a range and recording the relative reconstruction error and the orthogonality measure $\\|I - V_k^\\top V_k\\|_F$ per iteration; if even at large sketch sizes the randomized solution error does not approach the deterministic LSQR error or the basis orthogonality degrades linearly with iteration count, the claimed accuracy equivalence fails.","tokens_in":1702,"feed_emoji":"🎲","tokens_out":6307,"duration_ms":87322,"temperature":0.7,"pith_summary":"This paper develops randomized Krylov subspace methods for large-scale linear inverse problems with rectangular forward operators. Instead of computing exact inner products during orthogonalization, the methods estimate them through low-dimensional sketches of the vectors, which can reduce the cost of building a basis. The centerpiece is a randomized Golub-Kahan algorithm that produces two partial factorizations of the matrix, and from them three solvers - rLSQR, rCGLS, and rLSMR - that approximate least-squares solutions. The paper also builds hybrid versions that combine these projections with Tikhonov regularization and select the regularization parameter adaptively. Numerical tests on image deblurring and seismic tomography indicate that, with sufficiently large sketch sizes, the randomized solvers track their deterministic counterparts, while very small sketch sizes cause erratic behavior.","feed_headline":"Randomized Golub-Kahan solves rectangular inverse problems","feed_subtitle":"Sketched inner products cut orthogonalization cost while hybrid variants keep Tikhonov parameter choice adaptive.","key_machinery":"The central object is the randomized Golub-Kahan (rGK) factorization, built by the randomized Gram-Schmidt process with two sketching matrices, $\\Theta^{(n)}$ for vectors in $\\mathbb{R}^n$ and $\\Theta^{(m)}$ for vectors in $\\mathbb{R}^m$. At iteration $k$ it maintains stacks of sketched basis vectors and computes coefficients by sketched inner products, yielding the two factorizations in (2.14). These factorizations carry the argument because each solver is obtained by imposing a different condition on the residual or normal-equations residual and then reading off the projected problem from $M_k$ and $T_{k+1}$; the hybrid variants regularize those projected problems with Tikhonov terms and evaluate parameter-choice functionals cheaply.","core_discovery":"The central claim is that the randomized Golub-Kahan algorithm, Algorithm 2.2, gives exact-arithmetic factorizations $AV_k = U_{k+1}M_k$ and $A^\\top U_{k+1} = V_{k+1}T_{k+1}$ for a rectangular matrix $A$, where the columns of $V$ and $U$ are orthogonal with respect to the sketched inner products and $M_k$, $T_{k+1}$ are upper Hessenberg and upper triangular. From these factorizations the paper derives new solvers rLSQR, rCGLS, and rLSMR that generalize the corresponding deterministic Krylov methods: rLSQR minimizes the sketched residual norm over the Krylov subspace, rLSMR minimizes the sketched normal-equations residual, and rCGLS enforces a different orthogonality condition that has no short recurrence. The same framework yields hybrid Tikhonov-regularized variants in which the regularization parameter can be chosen at each iteration by the discrepancy principle or weighted GCV, because the approximation subspace does not depend on that parameter.","pith_inferences":["A natural testable extension, not pursued in the paper, is an adaptive scheme that enlarges the sketch dimensions when the sketched basis vectors lose orthogonality, rather than fixing the sketch size a priori.","Because the paper leaves the theoretical properties of the low-rank approximations from rGK open, establishing those properties could open the same factorizations to matrix-function estimation or uncertainty quantification in Bayesian inverse problems.","The paper notes that storage grows linearly with iterations because short recurrences are lost; combining rGK with the recycling strategies already mentioned for hybrid methods could keep memory bounded in long runs.","The observed sensitivity to sketch size suggests that a practical implementation should monitor the sketched orthogonality and terminate before semiconvergence deterioration, an operational rule the paper does not fully specify."],"forward_implications":["Rectangular inverse problems, including underdetermined and overdetermined ones, can now be attacked with randomized Krylov methods, not only square systems as with rGMRES.","When the sketch dimensions are chosen large enough, rLSQR, rCGLS, and rLSMR reproduce the semiconvergence behavior of standard LSQR, CGLS, and LSMR, so early termination still acts as regularization.","Hybrid randomized methods inherit the key property of deterministic hybrids: the approximation subspace is independent of the regularization parameter, so $\\lambda$ can be tuned at each iteration at negligible extra cost.","For problems needing only few iterations, such as high noise or severe ill-conditioning, smaller sketch sizes suffice, widening the regime of computational savings.","The rGK projected matrices approximate the dominant singular values of $A$, with rCGLS and rLSMR giving better approximations than rLSQR, which affects both iterative and hybrid regularization behavior."],"supporting_citations":[{"why":"Supplies the randomized Gram-Schmidt process and the rArnoldi/rGMRES framework that the new rGK approach extends to rectangular matrices.","marker":"[4]"},{"why":"Provides the hybrid projection method framework and the regularization parameter selection strategies used for the hybrid randomized solvers.","marker":"[10]"},{"why":"Provides the test problem package used to generate the image deblurring and seismic tomography examples.","marker":"[16]"},{"why":"Supplies the semiconvergence and iterative regularization viewpoint that motivates early termination in the randomized solvers.","marker":"[20]"},{"why":"Provides the randomized numerical linear algebra foundations, including oblivious subspace embedding background for the sketching matrices.","marker":"[25]"},{"why":"Supplies the general projection-method framework used to characterize the optimality properties of rLSQR and rLSMR.","marker":"[32]"},{"why":"Presents a related randomized and inner-product free Krylov approach for inverse problems that the paper builds on and contrasts with its own method.","marker":"[24]"},{"why":"Documents numerical instability of sketch-and-solve approaches for ill-conditioned problems, motivating the iterative randomized Krylov route taken here.","marker":"[27]"}],"fun_headline_variants":["Randomized Golub-Kahan handles rectangular inverse problems","Sketch-and-solve Krylov for rectangular inverse problems","Randomized Golub-Kahan extends Krylov to rectangular matrices","Randomized Krylov adapts Tikhonov on the fly for inverse problems","rLSQR and rLSMR: randomized Krylov for rectangular inverse problems"],"cache_read_input_tokens":23936,"weakest_assumption_plain":"The load-bearing premise is that the sketching matrices accurately preserve inner products over the Krylov subspaces generated during iteration; if the sketch size is too small, the basis vectors are not sufficiently orthogonal and the iterations become erratic.","fun_headline_variants_meta":{"raw":{"variants":["Randomized Golub-Kahan handles rectangular inverse problems","Sketch-and-solve Krylov for rectangular inverse problems","Randomized Golub-Kahan extends Krylov to rectangular matrices","Randomized Krylov adapts Tikhonov on the fly for inverse problems","rLSQR and rLSMR: randomized Krylov for rectangular inverse problems"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001051,"raw_usage":{"total_tokens":4396,"prompt_tokens":907,"completion_tokens":3489,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":523,"completion_tokens_details":{"reasoning_tokens":3400}},"tokens_in":523,"tokens_out":3489,"duration_ms":21639,"temperature":1.0,"reasoning_tokens":3400,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T16:49:06.031672+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run Algorithm 2.2 on a rectangular inverse problem, such as the seismic tomography example with $A$ of size $131072\\times 65536$, while increasing the sketch dimensions $\\ell_m/\\ell_n$ across a range and recording the relative reconstruction error and the orthogonality measure $\\|I - V_k^\\top V_k\\|_F$ per iteration; if even at large sketch sizes the randomized solution error does not approach the deterministic LSQR error or the basis orthogonality degrades linearly with iteration count, the claimed accuracy equivalence fails.","supporting_citations":[{"cited_title":"Chung and S","cited_arxiv_id":null,"evidence_quote":"Provides the hybrid projection method framework and the regularization parameter selection strategies used for the hybrid randomized solvers."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the semiconvergence and iterative regularization viewpoint that motivates early termination in the randomized solvers."},{"cited_title":"Martinsson and J","cited_arxiv_id":null,"evidence_quote":"Provides the randomized numerical linear algebra foundations, including oblivious subspace embedding background for the sketching matrices."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Presents a related randomized and inner-product free Krylov approach for inverse problems that the paper builds on and contrasts with its own method."},{"cited_title":"Meier, Y","cited_arxiv_id":null,"evidence_quote":"Documents numerical instability of sketch-and-solve approaches for ill-conditioned problems, motivating the iterative randomized Krylov route taken here."}],"review_version":1}