Pith. sign in

REVIEW 4 major objections 5 minor 28 references

Inexact Generalized Golub-Kahan Methods for Large-Scale Bayesian Inverse Problems

T0 review · 4 major / 5 minor · reviewed 2026-08-12 · deepseek-v4-flash

Pith's one-line read An inexact generalized Golub-Kahan solver matches exact-method accuracy for Bayesian inverse problems with uncertain forward models.

desk verdict A plausible but underdeveloped combination of inexact and generalized Golub-Kahan, with fixable pseudocode bugs and a missing perturbation analysis; worth a round of peer review if the authors can tighten it. read the letter →

arxiv 2411.14409 v1 pith:GCNZOGWU submitted 2024-11-21 math.NA cs.NA

classification math.NAcs.NA MSC 65F2265F1065K1015A29
keywords BayesianinverseproblemsgeneralizedTikhonovregularizationinexactGolub-KahandecompositionhybridprojectionmethodsKrylovsubspaceMaterncovariancekernelcomputedtomographyreconstructionparameterselection
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Bayesian inverse problems become hard in two realistic ways: the forward operator may only be known approximately, and the prior covariance matrix may be too large or irregular to factor. This paper develops an inexact generalized Golub-Kahan (igenGK) decomposition that accepts matrix-vector products with $A$ and $A^T$ that carry errors, and it uses the covariance matrix $Q$ only through matrix-vector products, never through its square root or inverse. On top of this decomposition the paper builds a hybrid iterative solver that projects the generalized Tikhonov problem into a small least-squares problem and chooses the regularization parameter automatically at each iteration. The central claim, supported by computed-tomography experiments, is that this inexact hybrid method reaches reconstruction accuracy comparable to methods that assume the exact forward operator.

What carries the argument

The carrying object is the inexact generalized Golub-Kahan (igenGK) decomposition: an iterative two-sided orthogonalization that combines the inexact Golub-Kahan idea with the generalized Golub-Kahan idea. At each step it multiplies by $(A+E_k)Q$ and by $(A+F_{k+1})^T R^{-1}$, orthogonalizes against the existing basis with respect to the $R^{-1}$ and $Q$ inner products, and records the coefficients in an upper Hessenberg matrix $M_k$ and a lower triangular matrix $L_{k+1}$. These two matrices convert the original large problem into a small projected Tikhonov problem, so the algorithm can be driven by standard parameter-choice rules on the projected problem.

What would settle it

Take a small matrix $A$, a covariance $Q$, and a known true solution; run the proposed igenGK hybrid method with controlled deterministic error matrices $E_k$ and $F_k$ of increasing norm and compare its reconstruction to the exact generalized Golub-Kahan method on the same problem. If a moderate error level (say $\|E_k\|/\|A\|$ around $10^{-2}$) yields relative reconstruction error substantially worse than the exact counterpart, or if a non-random error pattern of the same size causes divergence where random-error experiments stayed stable, the paper's central claim is falsified.

Watch

Extended reading notes

Core claim

After the change of variables $x = Q^{-1}(s - \mu)$, the MAP estimate becomes a generalized Tikhonov problem in which the forward operator acts as $AQ$. The paper's central object is a decomposition that, at iteration $k$, produces vectors $U_{k+1}$, $V_k$ and matrices $M_k$, $L_{k+1}$ satisfying $(A+E_k)QV_k = U_{k+1}M_k$ and $(A+F_{k+1})^T R^{-1}U_{k+1} = V_{k+1}L_{k+1}^T$, with $U$ and $V$ orthonormal in the $R^{-1}$ and $Q$ inner products respectively. Because the projected solution is found from the small problem $\min_y \|M_k y - \beta e_1\|_2^2 + \lambda^2 \|y\|_2^2$, the method never needs to form or factor $Q$, and it absorbs per-iteration errors $E_k$, $F_k$ in the forward-model products. The paper reports that the orthogonality conditions hold to machine precision and that the resulting reconstruction errors track the exact generalized method closely, including in a CT experiment where the inexactness comes from uncertain projection angles rather than synthetic random noise.

Load-bearing premise

The method works only if the errors in the approximate matrix-vector products are small enough that the reduced problem built from those approximate products still represents the original reconstruction problem; the paper gives no explicit bound on this, so this fidelity is the load-bearing assumption.

Editorial extensions

If this is right

  • The method solves the generalized Tikhonov MAP problem when the forward operator is available only through noisy matrix-vector products, so no exact $A$ or factorization of $Q$ is required.
  • In the computed-tomography experiments, the inexact generalized hybrid method's relative reconstruction error stays close to the exact genGK method at the same iterations.
  • Using the discrepancy principle for the regularization parameter gives results similar to using the optimal parameter, without needing the true solution.
  • The method remains stable when inexactness comes from realistic uncertain projection angles, provided the initial angle error is not too large.
  • Because $Q$ is accessed only via matrix-vector products, smooth stationary covariance kernels whose products can be accelerated (for example by Fourier methods) can be used directly.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • A natural next step, not taken in the paper, is to prove a perturbation bound showing that the noisy projected problem stays close to the exact one when $\|E_k\|$ and $\|F_k\|$ are small; the paper's evidence is entirely empirical.
  • The decreasing-inexactness experiment suggests a self-correcting loop: use the current reconstruction to improve estimates of uncertain parameters such as projection angles, then rerun the inexact solve with smaller errors.
  • The same decomposition should apply to other operator-uncertainty problems such as blind deblurring or super-resolution, where the error is structured model mismatch rather than Gaussian random noise.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 5 minor

Summary. The paper proposes an inexact generalized Golub-Kahan (igenGK) decomposition for solving large-scale Bayesian inverse problems when the forward operator is available only through approximate matrix-vector products. The method combines the inexact Golub-Kahan approach of Gazzola and Landman with the generalized Golub-Kahan bidiagonalization of Chung and Saibaba, and it avoids explicit factorizations of the prior covariance matrix Q by using only MVPs. A hybrid scheme with discrepancy principle and WGCV is developed for automatic regularization parameter selection. Numerical experiments on a CT reconstruction problem compare the method against standard, generalized, and inexact counterparts, and a separate experiment considers inexactness arising from inaccurate projection angles.

Significance. The paper addresses a relevant and timely problem: Bayesian inverse problems with an uncertain forward model. The algebraic derivation of the inexact generalized decomposition (Section 4.1) is straightforward but sound, and the ability to work with Q only through matrix-vector products is a practical strength. The numerical results on a synthetic CT problem, including a realistic angle-inexactness scenario, suggest that the method can be effective. However, the central claim—that the inexact generalized hybrid method achieves accuracy comparable to its exact counterpart—is supported only by a single numerical example, and the paper provides no perturbation analysis to establish when the noisy projected problem is a faithful surrogate for the true Tikhonov problem. There is also a structural issue in the core algorithm and a mismatch in the numerical verification. The work is an incremental combination of existing techniques, but it could be a useful contribution if these gaps are addressed.

major comments (4)
  1. [Section 4.1, Algorithm 4.1] The k=1 branch of Algorithm 4.1 is incomplete: it computes v_1 (line 3) but never computes u_2, the scalar m_{2,1}, or the updated entry in M_2, yet the return statement outputs U_{k+1} = U_2 and M_{k+1}. Consequently, the algorithm as written does not produce a well-defined basis for the first iteration. Please provide the missing initialization or restructure the loop so that the update for u_{k+1} is performed for k=1 as well.
  2. [Section 4.2, Eq. (4.5)] The projected problem (4.5) is derived using the inexact relation (A+E_k)QV_k = U_{k+1}M_k, but no bound is given relating the solution y_k of (4.5) to the exact solution x_λ of (2.2) in terms of ||E_k||, ||F_k||, and λ. The paper's main claim (end of Section 1) that the inexact hybrid method achieves accuracy comparable to its exact counterpart is therefore unsupported beyond the single tested example, and the method's reliability for larger inexactness levels is unclear. Please add a perturbation analysis or at least a precise discussion of the conditions on E_k and F_k under which the surrogate problem is a faithful approximation.
  3. [Section 5.1, Table 1] The verification in Table 1 checks the identity A^T U_k = V_k L_k^T, but the decomposition derived in (4.2) is (A+F_{k+1})^T R^{-1} U_{k+1} = V_{k+1} L_{k+1}^T. The table tests a different relation: it omits the R^{-1} weighting and uses the exact transpose A^T rather than the inexact operator. As a result, the numerical verification does not actually validate the relation that the projected problem (4.5) relies on. Please verify the correct identity (4.2), using the inexact transpose and the appropriate inner products.
  4. [Section 4.3, Eq. (4.8)] The residual in the discrepancy principle is defined as r_k(λ) = A s_k(λ) - d = M_k y_k(λ) - β_1 e_1. This equality is only valid if E_k = 0, because the second expression is the residual of the perturbed projected problem. Consequently, the DP stopping rule uses a residual that is biased by the model error, and the paper does not analyze how this bias affects the chosen regularization parameter. Please derive the relation between the true residual norm and the projected residual norm under the inexactness model, and discuss the implications for the DP criterion.
minor comments (5)
  1. [Throughout] There are several typos: 'cooresponds' in the Figure 2 caption, 'correponds' in Section 5.1, and 'Matérn' is spelled inconsistently. Please proofread carefully.
  2. [Section 3.3, Eq. (3.9)] The recurrence for v_{k+1} uses α_{k+1} twice, which is inconsistent with the standard genGK relation and with (3.10). Please check the coefficients and correct the equation.
  3. [Section 4.1] The notation for E_k and F_k is confusing: E_k is used both for the per-iteration error matrix in (3.5) and for the accumulated error defined after (4.2). Please introduce distinct notation for these two objects.
  4. [Section 5.1] It is not clear how iGK is applied to the generalized problem (2.2); if it requires a transformation or preconditioning, please specify the exact implementation to allow reproducibility.
  5. [Section 5.3, Figure 4] The WGCV results are reported to be worse than DP and optimal, but the paper does not provide the chosen weight parameter or any tuning strategy; a brief discussion would help the reader assess the practical usefulness of WGCV in this setting.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the proposed igenGK method is an algorithmic construction whose projected problem (4.5) is an algebraic reformulation of the noisy objective (4.3), and the central accuracy claim is an empirical comparison, not a prediction forced by construction.

full rationale

The derivation of the projected problem in Section 4.2 is not circular: equation (4.5) follows exactly from (4.3)-(4.4) using the defining relation (A+E_k)QV_k = U_{k+1}M_k and the R^{-1}-orthonormality of U_{k+1}, so the reduced objective is equivalent to the noisy objective by construction, but that is a valid reformulation rather than an input-output circle. The paper's main claim, that igenHyBR achieves accuracy comparable to its exact counterpart, is supported by numerical experiments (Figures 2-6) against independently implemented GK/iGK/genGK methods and against the true solution; no fitted parameter is relabeled as a prediction. Self-citations (e.g., Chung & Saibaba [8], Chung & Gazzola [3], Chung et al. [6]) provide prior algorithms and parameter-selection tools, but the central claim does not rest solely on those citations. A legitimate weakness is the absence of a perturbation bound controlling ||E_k|| and ||F_k|| or linking the noisy projected solution to the exact Tikhonov solution; this is a completeness/correctness risk, not a circularity. Therefore no circular step is present, and the score is 0.

Assumptions & free parameters 7 free parameters · 6 assumptions · 0 invented entities

The method introduces no new physical or mathematical entities; it combines existing algorithmic ingredients (inexact Golub-Kahan and generalized Golub-Kahan) with a new norm-compatible projection. The listed free parameters are experimental choices for the CT test problem, not fitted constants in a derivation. The method's validity depends on the standard Bayesian linear inverse problem assumptions and on the availability of fast Q MVPs.

free parameters (7)
  • beta (MVP inexactness level) = 1e-2, 1e-4, 1e-6 in experiments
    Hand-chosen Gaussian noise level for E_k and F_k to simulate inexact MVPs; not fitted to data.
  • nu (Matern kernel smoothness) = 1.5
    Hand-chosen hyperparameter of the prior covariance kernel.
  • ell (Matern correlation length) = 0.01
    Hand-chosen hyperparameter of the prior covariance kernel.
  • sigma (noise level) = set so ||epsilon||/||d_true|| = 0.04
    Noise variance in R = sigma I, chosen to produce 4% relative noise.
  • alpha_k (angle error schedule) = alpha_1 = 1e-1 or 1, decreasing logarithmically to 1e-6
    Hand-chosen schedule to simulate decreasing projection-angle errors in Section 5.4.
  • nu_DP (discrepancy principle constant) = not specified
    User-chosen constant in Eq. (4.8); value not reported.
  • omega (WGCV weight) = not specified
    User-chosen or estimated weight in Eq. (4.9); value not reported.
assumptions (6)
  • domain assumption The inverse problem is linear with additive Gaussian noise: d = A s + epsilon, epsilon ~ N(0,R).
    Basis of the Bayesian formulation in Section 1.
  • domain assumption The prior on s is Gaussian: s ~ N(mu, lambda^-2 Q).
    Used to derive the MAP estimate and generalized Tikhonov problem (1.3).
  • domain assumption A and Q are only accessible via matrix-vector products; Q is SPD and from the Matern family with fast (FFT/Toeplitz) MVPs.
    Core computational assumption in Sections 2 and 4.
  • domain assumption The inexact MVPs are of the form (A+E_k) and (A+F_k)^T with additive errors.
    The error model used throughout Sections 3.2 and 4.1.
  • standard math Standard Golub-Kahan bidiagonalization properties and Krylov subspace theory hold.
    The method is built on known results from [13, 21, 22].
  • standard math Bayes' theorem and the equivalence between MAP estimation and generalized Tikhonov regularization.
    Foundation of the problem formulation.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Inexact Generalized Golub-Kahan Methods for Large-Scale Bayesian Inverse Problems." pith.science (2026). https://pith.science/paper/GCNZOGWU

@misc{pith2026241114409,
  author       = {Pith},
  title        = {Pith review of: Inexact Generalized Golub-Kahan Methods for Large-Scale Bayesian Inverse Problems},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/GCNZOGWU}},
  note         = {Machine review of arXiv:2411.14409}
}
read the original abstract

Solving large-scale Bayesian inverse problems presents significant challenges, particularly when the exact (discretized) forward operator is unavailable. These challenges often arise in image processing tasks due to unknown defects in the forward process that may result in varying degrees of inexactness in the forward model. Moreover, for many large-scale problems, computing the square root or inverse of the prior covariance matrix is infeasible such as when the covariance kernel is defined on irregular grids or is accessible only through matrix-vector products. This paper introduces an efficient approach by developing an inexact generalized Golub-Kahan decomposition that can incorporate varying degrees of inexactness in the forward model to solve large-scale generalized Tikhonov regularized problems. Further, a hybrid iterative projection scheme is developed to automatically select Tikhonov regularization parameters. Numerical experiments on simulated tomography reconstructions demonstrate the stability and effectiveness of this novel hybrid approach.

Figures

Figures reproduced from arXiv: 2411.14409 by the authors.

Figure 1
Figure 1. Set up of parallel-beam X-ray CT where the detector collects data from a finite [PITH_FULL_IMAGE:figures/full_fig_p010_1.png] view at source ↗
Figure 2
Figure 2. Both (a) and (b) present comparisons of standard, generalized, inexact, and inexact [PITH_FULL_IMAGE:figures/full_fig_p012_2.png] view at source ↗
Figure 3
Figure 3. Comparisons of standard, generalized, inexact, and inexact generalized hybrid ap [PITH_FULL_IMAGE:figures/full_fig_p013_3.png] view at source ↗
Figures from the paper (3 more)
Figure 4
Figure 4. Figure 4: Comparison of optimal, DP, and WGCV methods in choosing the Tikhonov regu [PITH_FULL_IMAGE:figures/full_fig_p014_4.png]
Figure 5
Figure 5. Figure 5: Comparisons of reconstruction results for sinograms acquired at inexact angles, [PITH_FULL_IMAGE:figures/full_fig_p016_5.png]
Figure 6
Figure 6. Figure 6: Comparisons of reconstruction results for sinograms acquired at inexact angles, [PITH_FULL_IMAGE:figures/full_fig_p017_6.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

28 extracted references · 24 canonical work pages

  1. [1]

    Arioli, Generalized Golub–Kahan bidiagonalization and stopping criteria , SIAM Journal on Matrix Analysis and Applications, 34 (2013), pp

    M. Arioli, Generalized Golub–Kahan bidiagonalization and stopping criteria , SIAM Journal on Matrix Analysis and Applications, 34 (2013), pp. 571–592, https://doi.org/10.1137/120866543

  2. [2]

    Calvetti and E

    D. Calvetti and E. Somersalo , Bayesian scientific computing , Springer, 2023

  3. [3]

    Chung and S

    J. Chung and S. Gazzola, Computational methods for large-scale inverse problems: a survey on hybrid 17 BU projection methods, SIAM Review, 66 (2024), pp. 205–284

  4. [4]

    Chung, E

    J. Chung, E. Haber, and J. Nagy , Numerical methods for coupled super-resolution, Inverse Problems, 22 (2006), p. 1261

  5. [5]

    Chung and J

    J. Chung and J. G. Nagy , An efficient iterative approach for large-scale separable nonlinear inverse problems, SIAM Journal on Scientific Computing, 31 (2010), pp. 4654–4674

  6. [6]

    Chung, J

    J. Chung, J. G. Nagy, D. P. O’leary, et al. , A weighted GCV method for Lanczos hybrid regular- ization, Electronic Transactions on Numerical Analysis, 28 (2008), p. 2008

  7. [7]

    Chung and K

    J. Chung and K. Palmer , A hybrid LSMR algorithm for large-scale Tikhonov regularization , SIAM Journal on Scientific Computing, 37 (2015), pp. S562–S580

  8. [8]

    Chung and A

    J. Chung and A. K. Saibaba , Generalized hybrid iterative methods for large-scale Bayesian inverse problems, SIAM Journal on Scientific Computing, 39 (2017), pp. S24–S46, https://doi.org/10.1137/ 16M1081968

Show all 28 references
  1. [9]

    Elfving and P

    T. Elfving and P. C. Hansen, Unmatched projector/backprojector pairs: Perturbation and convergence analysis, SIAM Journal on Scientific Computing, 40 (2018), pp. A573–A591

  2. [10]

    Gazzola, P

    S. Gazzola, P. C. Hansen, and J. G. Nagy , IR Tools: a MATLAB package of iterative regularization methods and large-scale test problems , Numerical Algorithms, 81 (2019), pp. 773–811

  3. [11]

    Gazzola and M

    S. Gazzola and M. S. Landman , Regularization by inexact Krylov methods with applications to blind deblurring, SIAM Journal on Matrix Analysis and Applications, 42 (2021), pp. 1528–1552, https: //doi.org/10.1137/21M1402066

  4. [12]

    Gazzola and J

    S. Gazzola and J. G. Nagy , Generalized Arnoldi–Tikhonov method for sparse reconstruction , SIAM Journal on Scientific Computing, 36 (2014), pp. B225–B247

  5. [13]

    Golub and W

    G. Golub and W. Kahan, Calculating the singular values and pseudo-inverse of a matrix , Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis, 2 (1965), pp. 205–224

  6. [14]

    Golub and V

    G. Golub and V. Pereyra , Separable nonlinear least squares: the variable projection method and its applications, Inverse problems, 19 (2003), p. R1

  7. [15]

    Hanke and P

    M. Hanke and P. C. Hansen , Regularization methods for large-scale problems , Surv. Math. Ind, 3 (1993), pp. 253–315

  8. [16]

    P. C. Hansen , Discrete inverse problems: Insight and algorithms , SIAM, 2010

  9. [17]

    Ji and K

    H. Ji and K. W ang , Robust image deblurring with an inaccurate blur kernel , IEEE Transactions on Image processing, 21 (2011), pp. 1624–1634

  10. [18]

    M. E. Kilmer, P. C. Hansen, and M. I. Espanol , A projection-based approach to general-form Tikhonov regularization, SIAM Journal on Scientific Computing, 29 (2007), pp. 315–330

  11. [19]

    M. E. Kilmer and D. P. O’Leary, Choosing regularization parameters in iterative methods for ill-posed problems, SIAM Journal on matrix analysis and applications, 22 (2001), pp. 1204–1221

  12. [20]

    Nowak, S

    W. Nowak, S. Tenkleve, and O. A. Cirpka , Efficient computation of linearized cross-covariance and auto-covariance matrices of interdependent quantities , Mathematical Geology, 35 (2003), pp. 53–66

  13. [21]

    C. C. Paige and M. A. Saunders , Algorithm 583: LSQR: Sparse linear equations and least squares problems, ACM Transactions on Mathematical Software (TOMS), 8 (1982), pp. 195–209

  14. [22]

    C. C. Paige and M. A. Saunders , LSQR: An algorithm for sparse linear equations and sparse least squares, ACM Transactions on Mathematical Software (TOMS), 8 (1982), pp. 43–71

  15. [23]

    Reichel, F

    L. Reichel, F. Sgallari, and Q. Ye , Tikhonov regularization based on generalized Krylov subspace methods, Applied Numerical Mathematics, 62 (2012), pp. 1215–1228

  16. [24]

    R. A. Renaut, I. Hn ˇetynkov´a, and J. Mead , Regularization parameter estimation for large-scale Tikhonov regularization using a priori information , Computational statistics & data analysis, 54 (2010), pp. 3430–3445

  17. [25]

    Simoncini and D

    V. Simoncini and D. B. Szyld, Theory of inexact Krylov subspace methods and applications to scientific computing, SIAM Journal on Scientific Computing, 25 (2003), pp. 454–477

  18. [26]

    R. C. Smith , Uncertainty quantification: theory, implementation, and applications , SIAM, 2024

  19. [27]

    Uribe, J

    F. Uribe, J. M. Bardsley, Y. Dong, P. C. Hansen, and N. A. Riis, A hybrid Gibbs sampler for edge- preserving tomographic reconstruction with uncertain view angles, SIAM/ASA Journal on Uncertainty Quantification, 10 (2022), pp. 1293–1320

  20. [28]

    F. G. W aqar, S. Patel, and C. M. Simon, A tutorial on the Bayesian statistical approach to inverse problems, APL Machine Learning, 1 (2023), p. 041101, https://doi.org/10.1063/5.0154773, https: //doi.org/10.1063/5.0154773. 18

Pith tools

Reviewed August 12, 2026 · model on record in the stance chip above.