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 →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
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.
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
- 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.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
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)
- [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.
- [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.
- [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.
- [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)
- [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.
- [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.
- [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.
- [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.
- [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
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
free parameters (7)
- beta (MVP inexactness level) =
1e-2, 1e-4, 1e-6 in experiments
- nu (Matern kernel smoothness) =
1.5
- ell (Matern correlation length) =
0.01
- sigma (noise level) =
set so ||epsilon||/||d_true|| = 0.04
- alpha_k (angle error schedule) =
alpha_1 = 1e-1 or 1, decreasing logarithmically to 1e-6
- nu_DP (discrepancy principle constant) =
not specified
- omega (WGCV weight) =
not specified
assumptions (6)
- domain assumption The inverse problem is linear with additive Gaussian noise: d = A s + epsilon, epsilon ~ N(0,R).
- domain assumption The prior on s is Gaussian: s ~ N(mu, lambda^-2 Q).
- 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.
- domain assumption The inexact MVPs are of the form (A+E_k) and (A+F_k)^T with additive errors.
- standard math Standard Golub-Kahan bidiagonalization properties and Krylov subspace theory hold.
- standard math Bayes' theorem and the equivalence between MAP estimation and generalized Tikhonov regularization.
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 from the paper (3 more)
Reference graph
Works this paper leans on
-
[1]
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]
D. Calvetti and E. Somersalo , Bayesian scientific computing , Springer, 2023
work page 2023
-
[3]
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
work page 2024
- [4]
-
[5]
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
work page 2010
- [6]
-
[7]
J. Chung and K. Palmer , A hybrid LSMR algorithm for large-scale Tikhonov regularization , SIAM Journal on Scientific Computing, 37 (2015), pp. S562–S580
work page 2015
-
[8]
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
work page 2017
Show all 28 references
-
[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
2018
-
[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
2019
-
[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
2021 doi
-
[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
2014
-
[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
1965
-
[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
2003
-
[15]
Hanke and P
M. Hanke and P. C. Hansen , Regularization methods for large-scale problems , Surv. Math. Ind, 3 (1993), pp. 253–315
1993
-
[16]
P. C. Hansen , Discrete inverse problems: Insight and algorithms , SIAM, 2010
2010
-
[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
2011
-
[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
2007
-
[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
2001
-
[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
2003
-
[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
1982
-
[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
1982
-
[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
2012
-
[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
2010
-
[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
2003
-
[26]
R. C. Smith , Uncertainty quantification: theory, implementation, and applications , SIAM, 2024
2024
-
[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
2022
-
[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
2023 doi
Reviewed August 12, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.