pith. sign in

arxiv: 2502.03281 · v2 · pith:GJHRPG2Tnew · submitted 2025-02-05 · 🧮 math.NA · cs.NA

Efficient sampling approaches based on generalized Golub-Kahan methods for large-scale hierarchical Bayesian inverse problems

Elle Buser , Julianne Chung This is my paper

Pith reviewed 2026-05-23 04:03 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords hierarchical Bayesian inverse problemsMetropolis-Hastings samplingGolub-Kahan methodsuncertainty quantificationlarge-scale inverse problemsMCMC methodsseismic imagingphotoacoustic tomography
0
0 comments X

The pith

Generalized Golub-Kahan methods generate proposal distributions that enable efficient Metropolis-Hastings sampling inside Gibbs for hierarchical Bayesian inverse problems.

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

The paper develops computational techniques for sampling posterior distributions in large-scale linear inverse problems when both the unknown field and the hyperparameters controlling noise and prior variance must be inferred from data. Direct sampling is intractable because the conditional posterior for the field is high-dimensional and the square root and inverse of the prior covariance cannot be formed explicitly. The approach embeds Metropolis-Hastings independence sampling within a Gibbs loop, constructing the proposal from two variants of generalized Golub-Kahan bidiagonalization: a low-rank covariance approximation and a preconditioned Lanczos iteration. Numerical tests on seismic, photoacoustic, and atmospheric problems show that the resulting chains mix adequately while avoiding the prohibitive matrix operations.

Core claim

For hierarchical models with Gaussian priors and noise whose variances are treated as unknown hyperparameters, the conditional posterior for the unknown field remains Gaussian but high-dimensional. Generalized Golub-Kahan iterations produce low-rank or Lanczos-based approximations to this conditional covariance that serve as the covariance of an independence proposal inside Metropolis-Hastings; the acceptance probability is computed without ever forming the prior covariance square root or inverse. The resulting sampler therefore scales to problems where the number of unknowns is very large and direct covariance factorizations are infeasible.

What carries the argument

Generalized Golub-Kahan bidiagonalization applied to the forward operator and prior covariance to build low-rank or Lanczos approximations used as the covariance of the Metropolis-Hastings independence proposal for the conditional field posterior.

If this is right

  • Hyperparameters for noise variance and prior variance can be updated jointly with the field without forming full covariance matrices.
  • The method applies directly to dynamic inverse problems whose dimension grows with the number of time steps.
  • Two concrete proposal constructions (low-rank and preconditioned Lanczos) are supplied, each with its own cost-accuracy trade-off.
  • The same framework handles the three application classes shown: seismic imaging, dynamic photoacoustic tomography, and atmospheric inverse modeling.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same Golub-Kahan proposal construction could be inserted into other MCMC kernels, such as Hamiltonian Monte Carlo, whenever a Gaussian conditional appears.
  • Because the approach never requires an explicit prior covariance factorization, it may extend to priors whose covariance is only available through matrix-vector products.
  • If the forward operator changes slowly across iterations, the bidiagonalization factors could be warm-started, further reducing cost.

Load-bearing premise

The approximations produced by the generalized Golub-Kahan iterations remain accurate enough that the Metropolis-Hastings acceptance rate stays high enough for the chain to mix in reasonable time.

What would settle it

On a moderate-sized test problem where the exact conditional covariance can be formed directly, run both the proposed sampler and an exact Gaussian sampler from the same starting point and check whether the empirical distributions of the generated samples agree to within Monte-Carlo error.

read the original abstract

Uncertainty quantification for large-scale inverse problems remains a challenging task. For linear inverse problems with additive Gaussian noise and Gaussian priors, the posterior is Gaussian but sampling can be challenging, especially for problems with a very large number of unknown parameters (e.g., dynamic inverse problems) and for problems where computation of the square root and inverse of the prior covariance matrix are not feasible. Moreover, for hierarchical problems where several hyperparameters that define the prior and the noise model must be estimated from the data, the posterior distribution may no longer be Gaussian, even if the forward operator is linear. Performing large-scale uncertainty quantification for these hierarchical settings requires new computational techniques. In this work, we consider a hierarchical Bayesian framework where both the noise and prior variance are modeled as hyperparameters. Our approach uses Metropolis-Hastings independence sampling within Gibbs where the proposal distribution is based on generalized Golub-Kahan methods. We consider two proposal samplers, one that uses a low-rank approximation to the conditional covariance matrix and another that uses a preconditioned Lanczos method. Numerical examples from seismic imaging, dynamic photoacoustic tomography, and atmospheric inverse modeling demonstrate the effectiveness of the described approaches.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

0 major / 3 minor

Summary. The paper proposes two sampling methods based on generalized Golub-Kahan bidiagonalization to generate independence proposals for Metropolis-Hastings steps inside a Gibbs sampler for hierarchical Bayesian linear inverse problems. The approaches target conditional Gaussian posteriors (given hyperparameters) while avoiding explicit formation of the prior covariance square root or inverse, which is infeasible for large-scale problems. Effectiveness is illustrated through numerical experiments on seismic imaging, dynamic photoacoustic tomography, and atmospheric inverse modeling.

Significance. If the claimed efficiency and acceptance rates hold, the work supplies practical tools for uncertainty quantification in large-scale hierarchical inverse problems where standard Cholesky-based sampling is intractable. The use of existing Golub-Kahan infrastructure for proposal generation is a natural and potentially reusable extension of techniques already common in regularization and least-squares solvers.

minor comments (3)
  1. The abstract and introduction should explicitly state the precise form of the hierarchical model (e.g., which hyperparameters are sampled and the exact conditional distributions) to make the Gibbs structure immediately clear.
  2. Numerical results would benefit from reporting effective sample sizes or autocorrelation times in addition to acceptance rates to quantify mixing improvement over baseline methods.
  3. Clarify whether the low-rank and preconditioned-Lanczos proposals are used in separate experiments or compared head-to-head on the same problems.

Simulated Author's Rebuttal

0 responses · 0 unresolved

We thank the referee for the positive summary, recognition of the practical significance for large-scale hierarchical inverse problems, and recommendation of minor revision. No specific major comments were listed in the report.

Circularity Check

0 steps flagged

No significant circularity; algorithmic proposal is self-contained

full rationale

The paper presents two concrete algorithmic constructions (low-rank covariance approximation and preconditioned Lanczos) that extend standard Golub-Kahan bidiagonalization to generate independence proposals inside a Gibbs sampler for hierarchical Bayesian inverse problems. These constructions are described directly in terms of matrix-vector products with the forward operator and prior covariance without ever forming its square root or inverse; the acceptance rates and mixing are then checked on three independent numerical examples. No step reduces a claimed prediction or uniqueness result to a fitted quantity defined by the same parameters, nor does any load-bearing premise rest on a self-citation whose content is itself unverified. The derivation chain therefore remains externally falsifiable and does not collapse by construction.

Axiom & Free-Parameter Ledger

0 free parameters · 0 axioms · 0 invented entities

Based on abstract only, no specific free parameters, axioms, or invented entities can be identified from the provided information.

pith-pipeline@v0.9.0 · 5734 in / 1295 out tokens · 66329 ms · 2026-05-23T04:03:43.209373+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

24 extracted references · 24 canonical work pages

  1. [1]

    SIAM Review 66(2), 205–284 (2024)

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

    work page 2024

  2. [2]

    Geoscientific Model Devel- opment 15(14), 5547–5565 (2022) 27

    Cho, T., Chung, J., Miller, S.M., Saibaba, A.K.: Computationally efficient methods for large-scale atmospheric inverse modeling. Geoscientific Model Devel- opment 15(14), 5547–5565 (2022) 27

    work page 2022

  3. [3]

    Electronic Transactions on Numerical Analysis 58, 486–516 (2023)

    Pasha, M., Saibaba, A.K., Gazzola, S., Espa˜ nol, M.I., de Sturler, E.: A computa- tional framework for edge-preserving regularization in dynamic inverse problems. Electronic Transactions on Numerical Analysis 58, 486–516 (2023)

    work page 2023

  4. [5]

    Calvetti, D., Somersalo, E.: An Introduction to Bayesian Scientific Computing: Ten Lectures on Subjective Computing vol. 2. Springer, New York (2007)

    work page 2007

  5. [6]

    Bardsley, J.M.: Computational Uncertainty Quantification for Inverse Problems vol. 19. SIAM, Philadelphia (2018)

    work page 2018

  6. [7]

    SIAM Journal on Scientific Computing 35(6), 2494–2523 (2013)

    Bui-Thanh, T., Ghattas, O., Martin, J., Stadler, G.: A computational framework for infinite-dimensional Bayesian inverse problems part i: The linearized case, with application to global seismic inversion. SIAM Journal on Scientific Computing 35(6), 2494–2523 (2013)

    work page 2013

  7. [8]

    SIAM/ASA Jour- nal on Uncertainty Quantification 2(1), 203–222 (2014) https://doi.org/10.1137/ 120894877

    Bui-Thanh, T., Ghattas, O.: An analysis of infinite dimensional Bayesian inverse shape acoustic scattering and its numerical approximation. SIAM/ASA Jour- nal on Uncertainty Quantification 2(1), 203–222 (2014) https://doi.org/10.1137/ 120894877

    work page 2014

  8. [9]

    SIAM/ASA Journal on Uncertainty Quantification 6(3), 1076–1100 (2018)

    Brown, D.A., Saibaba, A., Vall´ elian, S.: Low-rank independence samplers in hierarchical Bayesian inverse problems. SIAM/ASA Journal on Uncertainty Quantification 6(3), 1076–1100 (2018)

    work page 2018

  9. [10]

    Numerical Linear Algebra with Applications 27(5), 2325 (2020)

    Saibaba, A.K., Chung, J., Petroske, K.: Efficient Krylov subspace methods for uncertainty quantification in large Bayesian linear inverse problems. Numerical Linear Algebra with Applications 27(5), 2325 (2020)

    work page 2020

  10. [11]

    Chapman & Hall/CRC Texts in Statisti- cal Science

    Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B.: Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statisti- cal Science. Taylor & Francis, Florida (2013). https://books.google.com/books? id=ZXL6AQAAQBAJ

    work page 2013

  11. [12]

    SIAM/ASA Journal on Uncertainty Quantification 7(3), 1105–1131 (2019)

    Saibaba, A.K., Bardsley, J., Brown, D.A., Alexanderian, A.: Efficient marginalization-based MCMC methods for hierarchical Bayesian inverse prob- lems. SIAM/ASA Journal on Uncertainty Quantification 7(3), 1105–1131 (2019)

    work page 2019

  12. [13]

    Springer, New York (2004)

    Robert, C.P., Casella, G.: Monte Carlo Statistical Methods, 2nd edn. Springer, New York (2004)

    work page 2004

  13. [14]

    Numerical Linear Algebra with Applications 27 (2020) https://doi.org/10.1002/nla.2325 28

    Saibaba, A., Chung, J., Petroske, K.: Efficient Krylov subspace methods for uncer- tainty quantification in large bayesian linear inverse problems. Numerical Linear Algebra with Applications 27 (2020) https://doi.org/10.1002/nla.2325 28

    work page doi:10.1002/nla.2325 2020

  14. [15]

    Inverse Problems 34(2), 024005 (2018)

    Chung, J., Saibaba, A.K., Brown, M., Westman, E.: Efficient generalized Golub– Kahan based methods for dynamic inverse problems. Inverse Problems 34(2), 024005 (2018)

    work page 2018

  15. [16]

    SIAM Journal on Matrix Analysis and Applications 34(2), 571–592 (2013)

    Arioli, M.: Generalized Golub–Kahan bidiagonalization and stopping criteria. SIAM Journal on Matrix Analysis and Applications 34(2), 571–592 (2013)

    work page 2013

  16. [17]

    SIAM Journal on Scientific Computing 39(5), 24–46 (2017)

    Chung, J., Saibaba, A.K.: Generalized hybrid iterative methods for large-scale Bayesian inverse problems. SIAM Journal on Scientific Computing 39(5), 24–46 (2017)

    work page 2017

  17. [18]

    Advances in Computational Mathematics 50(6), 1–33 (2024)

    Hall-Hooper, K.A., Saibaba, A.K., Chung, J., Miller, S.M.: Efficient iterative methods for hyperparameter estimation in large-scale linear inverse problems. Advances in Computational Mathematics 50(6), 1–33 (2024)

    work page 2024

  18. [19]

    IEEE Transactions on Information Theory 41(3), 613–627 (1995)

    Donoho, D.L.: De-noising by soft-thresholding. IEEE Transactions on Information Theory 41(3), 613–627 (1995)

    work page 1995

  19. [20]

    Numerical Algorithms 81(3), 773–811 (2019)

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

    work page 2019

  20. [21]

    Geoscientific Model Development 13(3), 1771–1785 (2020) https://doi.org/10.5194/gmd-13-1771-2020

    Miller, S.M., Saibaba, A.K., Trudeau, M.E., Mountain, M.E., Andrews, A.E.: Geostatistical inverse modeling with very large datasets: an example from the orbiting carbon observatory 2 (oco-2) satellite. Geoscientific Model Development 13(3), 1771–1785 (2020) https://doi.org/10.5194/gmd-13-1771-2020

    work page doi:10.5194/gmd-13-1771-2020 2020

  21. [22]

    Geoscientific Model Development 14(7), 4683–4696 (2021) https://doi.org/10.5194/gmd-14-4683-2021

    Liu, X., Weinbren, A.L., Chang, H., Tadi´ c, J.M., Mountain, M.E., Trudeau, M.E., Andrews, A.E., Chen, Z., Miller, S.M.: Data reduction for inverse modeling: an adaptive approach v1.0. Geoscientific Model Development 14(7), 4683–4696 (2021) https://doi.org/10.5194/gmd-14-4683-2021

    work page doi:10.5194/gmd-14-4683-2021 2021

  22. [23]

    Journal of Geophysical Research: Atmospheres 108(D16) (2003) https://doi.org/10.1029/2002JD003161 https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2002JD003161

    Lin, J.C., Gerbig, C., Wofsy, S.C., Andrews, A.E., Daube, B.C., Davis, K.J., Grainger, C.A.: A near-field tool for simulating the upstream influence of atmospheric observations: The stochastic time-inverted lagrangian transport (stilt) model. Journal of Geophysical Research: Atmospheres 108(D16) (2003) https://doi.org/10.1029/2002JD003161 https://agupubs....

    work page doi:10.1029/2002jd003161 2003

  23. [24]

    Meteorology and Atmospheric Physics 107, 51–64 (2010) https://doi.org/10.1007/s00703-010-0068-x

    Nehrkorn, T., Eluszkiewicz, J., Wofsy, S.C., Lin, J.C., Gerbig, C., Longo, M., Freitas, S.: Coupled weather research and forecasting–stochastic time-inverted lagrangian transport (wrf–stilt) model. Meteorology and Atmospheric Physics 107, 51–64 (2010) https://doi.org/10.1007/s00703-010-0068-x

    work page doi:10.1007/s00703-010-0068-x 2010

  24. [25]

    SIAM/ASA J

    Saibaba, A.K., Vallelian, S.: Low-rank independence samplers in hierarchical Bayesian inverse problems. SIAM/ASA J. Uncertain. Quantification6, 1076–1100 (2018) 29

    work page 2018