Computationally efficient methods for large-scale atmospheric inverse modeling

Atmospheric inverse modeling describes the process of estimating greenhouse gas fluxes or air pollution emissions at the Earth's surface using observations of these gases collected in the atmosphere. The launch of new satellites, the expansion of surface observation networks, and a desire for more detailed maps of surface fluxes have yielded numerous computational and statistical challenges for standard inverse modeling frameworks that were often originally designed with much smaller data sets in mind. In this article, we discuss computationally efficient methods for large-scale atmospheric inverse modeling and focus on addressing some of the main computational and practical challenges. We develop generalized hybrid projection methods, which are iterative methods for solving large-scale inverse problems, and specifically we focus on the case of estimating surface fluxes. These algorithms confer several advantages. They are efficient, in part because they converge quickly, they exploit efficient matrix–vector multiplications, and they do not require inversion of any matrices. These methods are also robust because they can accurately reconstruct surface fluxes, they are automatic since regularization or covariance matrix parameters and stopping criteria can be determined as part of the iterative algorithm, and they are flexible because they can be paired with many different types of atmospheric models. We demonstrate the benefits of generalized hybrid methods with a case study from NASA's Orbiting Carbon Observatory 2 (OCO-2) satellite. We then address the more challenging problem of solving the inverse model when the mean of the surface fluxes is not known a priori; we do so by reformulating the problem, thereby extending the applicability of hybrid projection methods to include hierarchical priors. We further show that by exploiting mathematical relations provided by the generalized hybrid method, we can efficiently calculate an approximate posterior variance, thereby providing uncertainty information.

1 Introduction

Numerous satellites and ground-based sensors observe greenhouse gas and air pollution mixing ratios in the atmosphere. A primary goal of atmospheric inverse modeling (AIM) is to estimate emissions or fluxes at the Earth's surface using these observations (Brasseur and Jacob, 2017; Enting, 2002; Michalak et al., 2004; Tarantola, 2005).

The number of greenhouse gas and air pollution measurements has greatly expanded in the past decade, enabling investigations of surface fluxes across larger regions, longer time periods, and/or at finer spatial and temporal detail. Carbon dioxide (CO₂) offers an illustrative example. The Greenhouse Gas Observing Satellite (GOSAT), the first satellite dedicated to monitoring CO₂ from space, launched in 2009 and collects $\sim \mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}$ high-quality observations globally each day (e.g., Nakajima et al., 2012). NASA's Orbiting Carbon Observatory 2 (OCO-2 satellite) launched in 2014 and collects ∼100 times more high-quality observations (Crisp, 2015; Eldering et al., 2017), and upcoming satellites like the Geostationary Carbon Observatory (GeoCarb) could collect up to $\sim \mathrm{1}\times {\mathrm{10}}^{\mathrm{7}}$ each day (though some of these observations will likely be unusable due to cloud cover) (Buis, 2018). These new observations are complemented by an expanding ground-based network of observations (NOAA Global Monitoring Laboratory, 2022) and expanded aircraft observations, including partnerships with several airlines to measure atmospheric CO₂ from regular commercial flights (Machida et al., 2008; Petzold et al., 2015). In addition, numerous cities are now monitored by dense urban networks of ground-based air sensors (e.g., Shusterman et al., 2016; Davis et al., 2017; Mitchell et al., 2018).

The growing spatial coverage and sheer number of atmospheric GHG observations make it increasingly possible to estimate GHG fluxes in greater spatial and temporal detail; it is key to developing inverse models that can represent all of the information in the observations by estimating emissions at sufficiently high spatiotemporal resolution. In addition, the increasing number of observation networks opens the possibility of rapid monitoring of the carbon cycle and monitoring of changes in anthropogenic emissions – from local to global scales. There is a concomitant need to report realistic uncertainty bounds on these emissions. Arguably, scientists require inverse modeling systems that can be used to estimate emissions and associated uncertainties quickly and at sufficient spatiotemporal resolution to facilitate this broad goal.

However, the computational challenges of achieving these goals are many. First, large-scale inverse models based on Bayesian statistics often require formulation of very large covariance matrices, calculation of matrix–matrix products with those covariance matrices, and/or solution of linear systems with those matrices. Second, existing inverse models often assume a Gaussian prior distribution for use with Bayes' theorem, where the prior mean vector and covariance matrix are required. Statistical approaches to estimating these covariance matrix parameters (e.g., restricted maximum likelihood estimation or Markov chain Monte Carlo methods) are often difficult to implement for extremely large inverse problems (Ganesan et al., 2014; Michalak et al., 2005), and a common approach is to populate the covariance matrices using expert knowledge. Third, fluxes often need to be estimated using iterative optimization algorithms for very large problems, and convergence of these algorithms can be slow (Miller et al., 2020). Fourth, calculating uncertainties in the estimated fluxes can be computationally prohibitive. The design of new methods to improve the computational feasibility of large atmospheric inverse problems has been the focus of numerous recent publications (see, e.g., Baker et al., 2006; Bousserez and Henze, 2018; Chatterjee and Michalak, 2013; Chatterjee et al., 2012; Chen et al., 2021a, b; Gourdji et al., 2012; Henze et al., 2007; Liu et al., 2021; Meirink et al., 2008; Miller et al., 2020, 2014; Yadav and Michalak, 2013; Zammit-Mangion et al., 2021).

Overview of features and contributions

The purpose of this study is to integrate several state-of-the-art computational and mathematical tools with AIM – tools that have been developed for and have had considerable success in other scientific fields (e.g., passive seismic tomography, medical imaging). Specifically, we investigate the use of generalized hybrid (genHyBR) projection methods for surface flux estimation and extend their use for inverse problems where the mean of the fluxes is not known a priori (sometimes referred to as geostatistical inverse modeling) (Chung and Saibaba, 2017; Saibaba et al., 2020). We address these challenges in our present work.

Building on prior work in Miller et al. (2020), we propose a unified computational framework for large-scale AIM with the following features.

1. We describe iterative genHyBR methods that are computationally efficient since they typically converge in a few iterations, are efficient in terms of storage, and work for very large satellite-based inverse problems. For example, we demonstrate the performance of genHyBR methods on two case studies previously considered in Miller et al. (2020). In the larger case study, we solve an inverse problem with 9×10⁶ unknown CO₂ fluxes and 1×10⁵ CO₂ observations.

2. We extend these methods to handle the case where the mean of the prior distribution is unknown, making genHyBR applicable to a broader range of inverse modeling applications that are common in the atmospheric science community (e.g., geostatistical inverse modeling).

3. Our approach is flexible in that it can be combined with any atmospheric transport model (e.g., either Lagrangian, particle-following models or the adjoint of an Eulerian model), and it can be used with a wide variety of covariance matrices for the unknown parameters and the noise.

4. Our framework also allows for efficiently estimating regularization parameters as part of the reconstruction, thus making it easier to objectively estimate the hyperparameters or covariance matrix parameters as part of the inverse model. We focus on the discrepancy principle (DP), which requires prior knowledge of an estimate of the noise, but provide alternate methods such as the unbiased predictive risk estimator and the generalized cross-validation, the latter of which does not require prior information regarding the noise level.

5. During the solution of the estimates, our solver stores information about the Krylov subspaces that can be used to estimate the posterior variance (at minimal computational cost), which gives insight into the uncertainty in the reconstructed solution. More precisely, evaluating uncertainties does not require additional model evaluations.

An overview of the paper is as follows. In Sect. 2, we describe the problem setup from a Bayesian perspective. In Sect. 3, we describe generalized hybrid methods for atmospheric inverse modeling. We show how to efficiently compute the maximum a posteriori (MAP) estimate and uncertainty estimates (e.g., posterior variance). The focus of this section is on the fixed mean case. We briefly mention an extension to the unknown mean case but defer most of the details to Appendix B. Numerical results are provided in Sect. 4, and conclusions can be found in Sect. 5.

2 Bayesian approach to inverse modeling

AIM will estimate greenhouse gas fluxes or air pollution emissions that match atmospheric observations, given an atmospheric transport model. It can be represented as an inverse problem of the form

$$\begin{array}{}\text{(1)}& \mathit{z}=\mathbf{H}\mathit{s}+\mathit{ϵ},\end{array}$$

where z∈ℝ^m is a vector of atmospheric observations, $\mathbf{H}\in {\mathbb{R}}^{m\times n}$ represents the forward atmospheric transport model, s∈ℝⁿ is a vector of the unknown surface fluxes or emissions, and ϵ∈ℝ^m represents noise or errors, including errors in the atmospheric observations (z) and in the atmospheric transport model (H). Note that s represents a vector containing spatial or spatiotemporal fluxes. Also, we assume $\mathit{ϵ}\sim \mathcal{N}(\mathbf{0},\mathbf{R})$, where $\mathbf{R}\in {\mathbb{R}}^{m\times m}$ is a positive definite matrix whose inverse and square root are inexpensive (e.g., a diagonal matrix with positive diagonal entries). The goal of the inverse problem is to estimate s given z and H. The inverse problem may be ill-posed or under-constrained by available observations. Therefore, it is common to include prior information to mitigate the ill-posedness, which is often referred to as variational regularization (Scherzer et al., 2009; Benning and Burger, 2018). We describe two different priors: fixed mean and unknown mean.

2.1 Fixed mean

A common approach in Bayesian inverse modeling is to model s as a Gaussian random variable with a fixed, known mean μ∈ℝⁿ and prior covariance matrix ${\mathit{\lambda }}^{-\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}\in {\mathbb{R}}^{n\times n}$. In many cases, this known mean (μ) is an emissions inventory, a bottom-up flux model, or a process-based model of CO₂ fluxes (Brasseur and Jacob, 2017). This approach is also known in the literature as Bayesian synthesis inversion. Using this framework, the prior distribution of s is given as follows:

$$\begin{array}{}\text{(2)}& \mathit{s}\sim \mathcal{N}(\mathit{\mu },{\mathit{\lambda }}^{-\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}).\end{array}$$

We assume that matrix Q_s is defined by a covariance kernel that describes the spatial and temporal variance and covariance in the prior distribution (Rasmussen and Williams, 2006). Furthermore, λ is a scaling parameter that is known a priori or has to be determined prior to the inversion process. The posterior distribution can be obtained by applying Bayes' theorem $\mathit{\pi }\left(\mathit{s}\right|\mathit{z})\propto \mathit{\pi }(\mathit{z}\left|\mathit{s}\right)\mathit{\pi }\left(\mathit{s}\right)$, which takes the form

$$\begin{array}{}\text{(3)}& \mathit{\pi }\left(\mathit{s}\right|\mathit{z})\propto \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert \mathbf{H}\mathit{s}-\mathit{z}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}-{\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \mathit{s}-\mathit{\mu }{\Vert }_{{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}}^{\mathrm{2}}\right),\end{array}$$

where $\Vert \mathit{x}{\Vert }_M^{\mathrm{2}}={\mathit{x}}^{\top }\mathbf{M}\mathit{x}$ for any symmetric positive definite matrix M and “∝” denotes a proportionality constant. The MAP estimate corresponding to this posterior distribution can be obtained by solving the optimization problem

$$\begin{array}{}\text{(4)}& {\mathit{s}}_{\mathrm{post}}:={\mathrm{argmin}}_{\mathit{s}\in {\mathbb{R}}^n}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert \mathbf{H}\mathit{s}-\mathit{z}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \mathit{s}-\mathit{\mu }{\Vert }_{{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}}^{\mathrm{2}}.\end{array}$$

Alternatively, it can be computed by solving the system of equations

$$({\mathbf{H}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathbf{H}+{\mathit{\lambda }}^{\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}){\mathit{s}}_{\mathrm{post}}={\mathbf{H}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathit{z}+{\mathit{\lambda }}^{\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathit{\mu }.$$

It is worth mentioning that the resulting posterior distribution is also Gaussian, with mean s_post and covariance Q_post, denoted as $\mathit{s}|\mathit{z}\sim \mathcal{N}({\mathit{s}}_{\mathrm{post}},{\mathbf{Q}}_{\mathrm{post}})$.

The reconstruction quality of Eq. (1) depends crucially on choosing appropriate covariance matrix parameters, or hyperparameters, that govern this prior Eq. (2) and the noise distribution of ϵ. In Sect. 3.1, we describe genHyBR methods for AIM where μ is fixed but λ is not known in advance.

In many applications, the prior mean μ may also not be known in advance and must be estimated as a part of the inversion process. Some inverse models (commonly referred to as geostatistical inverse models) directly assimilate environmental data or data on emitting activities directly into the inverse model, and the relationships between the surface fluxes and these data are rarely known a priori. In other cases, an emissions inventory or bottom-up flux model may be biased too high or too low. In these cases, Eq. (2) no longer holds, violating the statistical assumptions of the inverse model. One workaround is to scale the inventory or flux model as part of the inverse modeling process, which we now describe.

2.2 Unknown mean

In cases where the prior mean is unknown, we can represent the prior information in the form of the hierarchical model:

$$\begin{array}{}\text{(5)}& \mathit{s}|\mathit{\beta }\sim \mathcal{N}(\mathbf{X}\mathit{\beta },{\mathit{\lambda }}^{-\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}),\mathit{\beta }\sim \mathcal{N}({\mathit{\mu }}_{\mathit{\beta }},{\mathit{\lambda }}_{\mathit{\beta }}^{-\mathrm{2}}{\mathbf{Q}}_{\mathit{\beta }}).\end{array}$$

Here $\mathbf{X}\in {\mathbb{R}}^{n\times p}$ is a fixed matrix that includes covariates (e.g., environmental data or activity data) or a bottom-up inventory/flux model, ${\mathbf{Q}}_{\mathrm{s}}\in {\mathbb{R}}^{n\times n}$ is the prior covariance matrix, and λ is a scaling parameter. A set of unknown coefficients β∈ℝ^p scales the columns of X and is estimated as part of the inverse model. These coefficients are assumed to follow a Gaussian distribution with a given mean μ_β∈ℝ^p, covariance matrix ${\mathbf{Q}}_{\mathit{\beta }}\in {\mathbb{R}}^{p\times p}$, and scaling parameter λ_β.

Given the assumptions in Eqs. (1) and (5), from Bayes' theorem the posterior probability density function for the unknown mean case can be written as

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathit{\pi }(\mathit{s},\mathit{\beta }|\mathit{z})\propto & \mathit{\pi }\left(\mathit{z}\right|\mathit{s},\mathit{\beta }\left)\mathit{\pi }\right(\mathit{s}\left|\mathit{\beta }\right)\mathit{\pi }\left(\mathit{\beta }\right)\\ \propto & \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert \mathbf{H}\mathit{s}-\mathit{z}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}-{\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \mathit{s}-\mathbf{X}\mathit{\beta }{\Vert }_{{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}}^{\mathrm{2}}\right.\\ & \left.-{\displaystyle \frac{{\mathit{\lambda }}_{\mathit{\beta }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \mathit{\beta }-{\mathit{\mu }}_{\mathit{\beta }}{\Vert }_{{\mathbf{Q}}_{\mathit{\beta }}^{-\mathrm{1}}}^{\mathrm{2}}\right).\end{array}\end{array}$$

The MAP estimate can be written as the solution of the optimization problem

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\gamma }}_{\mathrm{post}}& ={\mathrm{argmin}}_{\mathit{\gamma }=[{\mathit{s}}^{\top },{\mathit{\beta }}^{\top }{]}^{\top }}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert \mathbf{H}\mathit{s}-\mathit{z}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}\\ & +{\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \mathit{s}-\mathbf{X}\mathit{\beta }{\Vert }_{{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\lambda }}_{\mathit{\beta }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \mathit{\beta }-{\mathit{\mu }}_{\mathit{\beta }}{\Vert }_{{\mathbf{Q}}_{\mathit{\beta }}^{-\mathrm{1}}}^{\mathrm{2}}.\end{array}\end{array}$$

The posterior distribution in Eq. (6) is Gaussian; therefore, the mean of the posterior distribution is also the MAP estimate, and the covariance is the inverse of the Hessian matrix of Eq. (7) and is given by

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\mathbf{\Gamma }}_{\mathrm{post}}=\\ & {\left[\begin{array}{cc}{\mathit{\lambda }}^{\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}+{\mathbf{H}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathbf{H}& -{\mathit{\lambda }}^{\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}\\ -{\mathit{\lambda }}^{\mathrm{2}}{\mathbf{X}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}& {\mathit{\lambda }}_{\mathit{\beta }}^{\mathrm{2}}{\mathbf{Q}}_{\mathit{\beta }}^{-\mathrm{1}}+{\mathit{\lambda }}^{\mathrm{2}}{\mathbf{X}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}\end{array}\right]}^{-\mathrm{1}}.\end{array}\end{array}$$

Therefore, the resulting posterior distribution is

$$\mathit{\gamma }\mid \mathit{z}\sim \mathcal{N}\left({\mathit{\gamma }}_{\mathrm{post}},{\mathbf{\Gamma }}_{\mathrm{post}}\right).$$

Here, λ and λ_β are scaling parameters that may not be known in advance, but we assume that λ_β=αλ with a constant α>0, where α is set in advance. We describe genHyBR methods for the unknown mean case in Sect. 3.3.

Note that previous works (Miller et al., 2020; Saibaba and Kitanidis, 2015) assume an improper prior for β (i.e., p(β)∝1), in which case a solution estimate can be obtained as $\widehat{\mathit{s}}=\mathbf{X}\widehat{\mathit{\beta }}+{\mathbf{Q}}_{\mathrm{s}}{\mathbf{H}}^{\top }\widehat{\mathit{\xi }},$ where

$$\begin{array}{}\text{(9)}& \left[\begin{array}{cc}{\mathbf{HQ}}_{\mathrm{s}}{\mathbf{H}}^{\top }+\mathbf{R}& \mathbf{HX}\\ (\mathbf{HX}{)}^{\top }& \mathbf{0}\end{array}\right]\left[\begin{array}{c}\widehat{\mathit{\xi }}\\ \widehat{\mathit{\beta }}\end{array}\right]=\left[\begin{array}{c}\mathit{z}\\ \mathbf{0}\end{array}\right].\end{array}$$

The system in Eq. (9) is often referred to as the dual-function form, and there are several equivalent formulations of these equations (Michalak et al., 2004). The size of the resulting system of equations is $(m+p)\times (m+p)$, where m is the number of measurements and p is the number of unknown parameters in β, so forming or inverting the matrix in Eq. (9) is unfeasible in many applications. The approach taken in Saibaba and Kitanidis (2012) and Miller et al. (2020) uses a matrix-free iterative method to solve Eq. (9); however, the number of required iterations can be very large, especially for problems with many measurements and even with the use of a preconditioner to speed convergence.

In this paper, we follow a different approach to handle the unknown mean case by using iterative hybrid approaches to a reformulated problem. Since these methods work directly on the least squares problem Eq. (7), the number of unknown parameters is n+p. However, the size of the linear system that defines the MAP estimate is independent of the number of observations, making it attractive for large data sets. Furthermore, our framework can handle a wide class of prior covariance operators, where the resulting prior covariance matrices are large and dense, and explicitly forming and factorizing these matrices is prohibitively expensive. These include, for example, prior covariance matrices that arise from nonseparable, spatiotemporal covariance kernels and parameterized kernels on non-uniform grids. Our approach only relies on forming matrix–vector products with the covariance matrices and is compatible with acceleration techniques using fast Fourier transform (FFT) or hierarchical matrix approaches. See Chung and Saibaba (2017) and Chung et al. (2018) for a detailed discussion. Thus, as we show in Sect. 4, our approach can incorporate various prior models and can scale to very large data sets.

3 Generalized hybrid projection methods for AIM

In this section, we describe generalized hybrid projection methods, dubbed genHyBR methods, for AIM. Hybrid projection methods were first developed in the 1980s as a way to combine iterative projection methods (e.g., Krylov subspace methods) and variational regularization methods (e.g., Tikhonov regularization) for solving very large inverse problems. These are iterative methods, where each iteration requires the expansion of the solution subspace, the estimation of the regularization parameter(s), and the solution of a projected, regularized problem. We point the interested reader to survey papers (Chung and Gazzola, 2021; Gazzola and Sabaté Landman, 2020). In Chung and Saibaba (2017), genHyBR methods were developed for computing Tikhonov-regularized solutions to problems where explicit computations of the square root and inverse of the prior covariance matrix are not feasible. This work enabled hybrid projection methods for more general regularization terms. The main benefits of genHyBR methods that make them ideal for large-scale AIM are efficiency, due to fast convergence to an accurate reconstruction of surface fluxes where efficient matrix–vector multiplications are exploited at each iteration, automatic estimation of parameters (e.g., hyperparameters and algorithmic parameters), and flexibility, because they can be paired with many different atmospheric models.

We describe genHyBR methods for both the fixed mean case (Sect. 3.1) and the unknown mean case (Sect. 3.3), with particular emphasis on the associated challenges for large data sets and subsequent uncertainty quantification (UQ) (Sect. 3.2). We provide a general overview of our approach, including the main components of genHyBR methods, in the flowchart in Fig. 1.

Figure 1 This flowchart provides a general overview of using genHyBR methods for AIM and subsequent UQ. Given input (corresponding to the observations and details of the problem setup), genHyBR is an iterative approach to approximate the MAP estimate. Each iteration of genHyBR consists of expanding the solution subspace (Sect. 3.1.1), projecting the problem, estimating a regularization parameter (Sect. 3.1.2), and solving a projected, regularized problem. After obtaining the MAP estimate, information computed from genHyBR can be used to efficiently estimate the posterior variance for UQ (Sect. 3.2).

3.1 Generalized hybrid methods with a fixed mean

To introduce genHyBR methods, we begin with the fixed mean case described in Sect. 2. If symmetric decompositions ${\mathbf{R}}^{-\mathrm{1}}={\mathbf{L}}_{\mathbf{R}}^{\top }{\mathbf{L}}_{\mathbf{R}}$ and ${\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}={\mathbf{L}}_{\mathrm{s}}^{\top }{\mathbf{L}}_{\mathrm{s}}$ are available, then optimization problem Eq. (4) can be rewritten in the standard least squares form

$${\mathit{s}}_{\mathrm{post}}={\mathrm{argmin}}_{\mathit{s}\in {\mathbb{R}}^n}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{∥\left[\begin{array}{c}{\mathbf{L}}_{\mathbf{R}}\mathbf{H}\\ \mathit{\lambda }{\mathbf{L}}_{\mathrm{s}}\end{array}\right]\mathit{s}-\left[\begin{array}{c}{\mathbf{L}}_{\mathbf{R}}\mathit{z}\\ \mathit{\lambda }{\mathbf{L}}_{\mathrm{s}}\mathit{\mu }\end{array}\right]∥}_{\mathrm{2}}^{\mathrm{2}}.$$

However, computing L_s can be computationally infeasible for large n, and λ may not be known a priori. This motivates us to use the following change in variables,

$$\begin{array}{}\text{(10)}& \mathit{x}\leftarrow {\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}(\mathit{s}-\mathit{\mu }),\mathit{b}\leftarrow \mathit{z}-\mathbf{H}\mathit{\mu },\end{array}$$

in which case the solution to problem Eq. (4) is given by ${\mathit{s}}_{\mathrm{post}}=\mathit{\mu }+{\mathbf{Q}}_{\mathrm{s}}\mathit{x}$, where x solves

$$\begin{array}{}\text{(11)}& \underset{\mathit{x}\in {\mathbb{R}}^n}{min}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert {\mathbf{HQ}}_{\mathrm{s}}\mathit{x}-\mathit{b}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \mathit{x}{\Vert }_{{\mathbf{Q}}_{\mathrm{s}}}^{\mathrm{2}}.\end{array}$$

Note that, with this reformulation, we avoid L_s, ${\mathbf{L}}_{\mathrm{s}}^{-\mathrm{1}}$, and ${\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}$ and only require matrix–vector products with Q_s. Furthermore, for iterative methods for Eq. (11), the matrix H does not need to be formed explicitly, as we only need access to matrix–vector and matrix–transpose–vector products.

There are two main ingredients in the genHyBR approach: (1) the generalized Golub–Kahan bidiagonalization approach for constructing a solution subspace and (2) regularization parameter estimation methods for computing a suitable regularization parameter in the projected space.

3.1.1 Generalized Golub–Kahan bidiagonalization

We now describe the generalized Golub–Kahan bidiagonalization approach that is the backbone of the genHyBR method. Given matrices H, R, and Q_s and vector b from Eq. (11), the basic idea is to generate a set of basis vectors contained in V_k for the Krylov subspace:

$${\mathcal{S}}_k\equiv \mathcal{R}\left({\mathbf{V}}_k\right)={\mathcal{K}}_k({\mathbf{H}}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathbf{HQ}}_{\mathrm{s}},{\mathbf{H}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathit{b}),$$

where ℛ(⋅) denotes the column space and the Krylov subspace is ${\mathcal{K}}_k(\mathbf{M},\mathit{f})=\text{Span}\mathit{\{}\mathit{f},\mathbf{M}\mathit{f},\mathrm{\dots },{\mathbf{M}}^{k-\mathrm{1}}\mathit{f}\mathit{\}}$. The generated basis vectors span a low-dimensional subspace that is rich in information about important directions; thus, solutions to the (smaller) projected problem often provide good approximations to the solution of the high-dimensional problem. With initializations ${\mathit{\delta }}_{\mathrm{1}}=\Vert \mathit{b}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}$, ${\mathit{u}}_{\mathrm{1}}=\mathit{b}/{\mathit{\delta }}_{\mathrm{1}}$ and ${\mathit{\gamma }}_{\mathrm{1}}{\mathit{v}}_{\mathrm{1}}={\mathbf{H}}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathit{u}}_{\mathrm{1}}$, the kth iteration of the generalized Golub–Kahan bidiagonalization procedure generates vectors u_k+1 and v_k+1 such that

$$\begin{array}{ll}{\displaystyle }{\mathit{\gamma }}_{k+\mathrm{1}}{\mathit{u}}_{k+\mathrm{1}}=& {\displaystyle }{\mathbf{HQ}}_{\mathrm{s}}{\mathit{v}}_k-{\mathit{\gamma }}_k{\mathit{u}}_k,\\ {\displaystyle }{\mathit{\delta }}_{k+\mathrm{1}}{\mathit{v}}_{k+\mathrm{1}}=& {\displaystyle }{\mathbf{H}}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathit{u}}_{k+\mathrm{1}}-{\mathit{\delta }}_{k+\mathrm{1}}{\mathit{v}}_k,\end{array}$$

where scalars ${\mathit{\gamma }}_k,{\mathit{\delta }}_k\ge \mathrm{0}$ are computed such that $\Vert {\mathit{u}}_k{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}=\Vert {\mathit{v}}_k{\Vert }_{{\mathbf{Q}}_{\mathrm{s}}}=\mathrm{1}$. At the end of k iterations, we have

$$\begin{array}{rl}{\mathbf{B}}_k& \equiv \left[\begin{array}{cccc}{\mathit{\gamma }}_{\mathrm{1}}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ {\mathit{\delta }}_{\mathrm{2}}& {\mathit{\gamma }}_{\mathrm{2}}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& {\mathit{\delta }}_{\mathrm{3}}& \mathrm{\ddots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{\ddots }& {\mathit{\gamma }}_k\\ \mathrm{0}& \mathrm{\cdots }& \mathrm{0}& {\mathit{\delta }}_{k+\mathrm{1}}\end{array}\right],\\ & {\mathbf{U}}_{k+\mathrm{1}}\equiv \left[\begin{array}{c}{\mathit{u}}_{\mathrm{1}},\mathrm{\dots },{\mathit{u}}_{k+\mathrm{1}}\end{array}\right],\text{ and }{\mathbf{V}}_k\equiv \left[\begin{array}{c}{\mathit{v}}_{\mathrm{1}},\mathrm{\dots },{\mathit{v}}_k\end{array}\right],\end{array}$$

where the following relations hold:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\mathbf{HQ}}_{\mathrm{s}}{\mathbf{V}}_k& ={\mathbf{U}}_{k+\mathrm{1}}{\mathbf{B}}_k\text{and}\\ & {\mathbf{H}}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathbf{U}}_{k+\mathrm{1}}={\mathbf{V}}_k{\mathbf{B}}_k^{\top }+{\mathit{\gamma }}_{k+\mathrm{1}}{\mathit{v}}_{k+\mathrm{1}}{\mathit{e}}_{k+\mathrm{1}}^{\top },\end{array}\end{array}$$

where e_j is the jth column of an identity matrix with the appropriate dimensions. Also, matrices U_k+1 and V_k satisfy the following orthogonality conditions:

$$\begin{array}{}\text{(13)}& {\mathbf{U}}_{k+\mathrm{1}}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathbf{U}}_{k+\mathrm{1}}={\mathbf{I}}_{k+\mathrm{1}}\text{ and }{\mathbf{V}}_k^{\top }{\mathbf{Q}}_{\mathrm{s}}{\mathbf{V}}_k={\mathbf{I}}_k,\end{array}$$

with ${\mathbf{U}}_{k+\mathrm{1}}{\mathit{\delta }}_{\mathrm{1}}{\mathit{e}}_{\mathrm{1}}=\mathit{b}$. Then, for a given λ>0, the solution to Eq. (11) is recovered by ${\mathit{x}}_{k,\mathit{\lambda }}={\mathbf{Q}}_{\mathrm{s}}{\mathbf{V}}_k{\mathit{y}}_{k,\mathit{\lambda }}$, where y_k,λ is the solution to the regularized, projected problem

$$\begin{array}{}\text{(14)}& \underset{\mathit{y}\in {\mathbb{R}}^k}{min}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert {\mathbf{B}}_k\mathit{y}-{\mathit{\delta }}_{\mathrm{1}}{\mathit{e}}_{\mathrm{1}}{\Vert }_{\mathrm{2}}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \mathit{y}{\Vert }_{\mathrm{2}}^{\mathrm{2}}.\end{array}$$

Note that Eq. (14) is a standard least squares problem with Tikhonov regularization, and since the coefficient matrix B_k is of size $(k+\mathrm{1})\times k$, the solution can be computed efficiently (Björck, 1996). Each iteration of the generalized Golub–Kahan bidiagonalization process requires one matrix–vector product with H and its adjoint (suppose we denote its cost by T_H), two matrix–vector products with Q_s (similarly, denoted $T_{{\mathbf{Q}}_{\mathrm{s}}}$), and additional 𝒪(m+n) floating point operations (flops). To compute the solution of the least squares problem Eq. (14), the cost is 𝒪(k³) flops, and the cost of forming x_k,λ is 𝒪(nk) flops. Thus, the overall cost of the algorithm is

$$T_{\text{genGK}}=\mathrm{2}k(T_{\mathbf{H}}+T_{{\mathbf{Q}}_{\mathrm{s}}})+\mathcal{O}\left(k\right(m+n\left)\right)\text{flops}.$$

In practice, the vectors {u_k} and {v_k} lose orthogonality in floating point arithmetic, so full or partial re-orthogonalization (Barlow, 2013) can be used to ensure orthogonality. This costs an additional 𝒪(k²(m+n)) flops. Thus far we have described an iterative method for approximating the MAP estimate, Eq. (4), where the kth iterate is given by ${\mathit{s}}_{k,\mathit{\lambda }}=\mathit{\mu }+{\mathbf{Q}}_{\mathrm{s}}{\mathit{x}}_{k,\mathit{\lambda }}$ for fixed regularization parameter λ.

3.1.2 Regularization parameter estimation methods

In this subsection, we highlight one of the main computational benefits of hybrid projection methods, which is the ability to estimate regularization parameters efficiently and adaptively while still ensuring robustness of the solution. For genHyBR approaches, we use the generalized Golub–Kahan bidiagonalization to generate a projection subspace and solve the projected problem Eq. (14) while simultaneously estimating the regularization parameter λ. Note that the regularization term in the projected system (14) is standard Tikhonov regularization, and a plethora of parameter estimation methods exist for Tikhonov regularization (Bardsley, 2018; Hansen, 2010). Here we focus on the discrepancy principle (DP) and point the interested reader to Appendix A for further details on other parameter estimation methods that can be incorporated within genHyBR methods for AIM.

The DP is a common approach for estimating a regularization parameter, where the main goal is to determine λ such that the residual norm for the regularized reconstruction matches a given estimate of the noise level in the observations. That is, the DP method selects the largest parameter value λ for which the reconstructed fluxes s_λ satisfy

$$\begin{array}{}\text{(15)}& {\mathcal{D}}_{\mathrm{full}}\left(\mathit{\lambda }\right)=\Vert \mathbf{H}{\mathit{s}}_{\mathit{\lambda }}-\mathit{z}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}\le \mathit{\tau }m,\end{array}$$

where τ≥1 is a user-defined parameter and m is the expected value of $\Vert \mathit{ϵ}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}$. Typical choices for safety factor τ are in the range $\mathrm{1}\le \mathit{\tau }\le \mathrm{2}.$

The DP has been used in AIM (Hase et al., 2017) as well as more generally in inverse problems (Groetsch, 1983; Hansen, 2010). For a given λ, evaluating 𝒟_full(λ) requires computation of s_λ and matrix–vector multiplication by H, which can get costly if many different values of λ are desired. However, a major distinction here is that by using an iterative hybrid formulation, we can exploit relationships in Eq. (12) for more efficient parameter selection. In particular, the residual norm can be simplified as

$$\begin{array}{}\text{(16)}& \Vert \mathbf{H}{\mathit{s}}_{k,\mathit{\lambda }}-\mathit{z}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}=\Vert {\mathbf{B}}_k{\mathit{y}}_{k,\mathit{\lambda }}-{\mathit{\delta }}_{\mathrm{1}}{\mathit{e}}_{\mathrm{1}}{\Vert }_{\mathrm{2}}^{\mathrm{2}}\equiv {\mathcal{D}}_{\mathrm{proj}}\left(\mathit{\lambda }\right).\end{array}$$

Thus, we let λ_k be the regularization parameter estimated for the projected problem at the kth iteration, such that ${\mathcal{D}}_{\mathrm{proj}}\left({\mathit{\lambda }}_k\right)\le \mathit{\tau }m.$ Then, as the number of iterations k increases, the estimated DP regularization parameter for the projected problem becomes a better approximation of the DP parameter for the original problem.

The advantage of this approach is two-fold: the regularization parameter is selected adaptively (i.e., each iteration can have a different regularization parameter), and the cost of parameter selection is cheap (𝒪(k³) flops) since we work with small matrices of size $(k+\mathrm{1})\times k$ and k is much smaller than m and n. Furthermore, there are various theoretical results that show that selecting the regularization parameter for the projected problem (i.e., project then regularize) is equivalent to first estimating the regularization parameter and then using an iterative projection method (i.e., regularize then project) (Chung and Gazzola, 2021).

3.2 Approximation to the posterior covariance matrices

In the Bayesian approach for fixed parameter λ, the posterior distribution Eq. (3) is Gaussian and, thus, is fully specified by the mean and covariance matrix. However, neither computing nor storing the covariance matrix is feasible, making further uncertainty estimation challenging. Instead, we follow the approach described in Chung et al. (2018) and Saibaba et al. (2020) for the fixed mean case, where an approximation to the posterior covariance matrix is obtained using the computed vectors generated during the generalized Golub–Kahan bidiagonalization process. An advantage of this approach is that, by storing partial information while computing the MAP estimate, we can approximately compute the uncertainty associated with the MAP estimate (e.g., posterior variance) with minimal additional cost and no further access to the forward and adjoint models. For the fixed mean case, we refer to this approach as genHyBRs with UQ and provide a summary in Algorithm 1.

Algorithm 1 AIM with fixed mean – genHyBRs with UQ.

Require:

Matrices H, R and Q, and vector b.

1:

{ Compute MAP estimate }

2:

initialize ${\mathit{u}}_{\mathrm{1}}=\mathit{b}/\Vert \mathit{b}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}$

3:

for $j=\mathrm{1},\mathrm{\dots },k$ do

4:

one iteration of generalized Golub–Kahan bidiagonalization to obtain ${\mathbf{B}}_j,{\mathbf{U}}_{j+\mathrm{1}},$ and V_j

5:

estimate regularization parameter λ and compute

${\mathit{x}}_{j,\mathit{\lambda }}={\mathbf{Q}}_{\mathrm{s}}{\mathbf{V}}_j{\mathit{y}}_{j,\mathit{\lambda }}$ where y_j,λ solves Eq. (14)

6:

end for

7:

compute the MAP estimate ${\mathit{s}}_{k,\mathit{\lambda }}=\mathit{\mu }+{\mathbf{Q}}_{\mathrm{s}}{\mathit{x}}_{k,\mathit{\lambda }}$

8:

{ Compute the approximation to the posterior variance }

9:

compute the eigendecomposition ${\mathbf{B}}_k^{\top }{\mathbf{B}}_k={\mathbf{W}}_k{\mathbf{\Theta }}_k{\mathbf{W}}_k^{\top }$

10:

compute Z_k=Q_sV_kW_k and

${\mathbf{\Delta }}_k={\mathit{\lambda }}^{-\mathrm{2}}\text{diag}(\frac{{\mathit{\theta }}_{\mathrm{1}}}{{\mathit{\theta }}_{\mathrm{1}}+{\mathit{\lambda }}^{\mathrm{2}}},\mathrm{\dots },\frac{{\mathit{\theta }}_k}{{\mathit{\theta }}_k+{\mathit{\lambda }}^{\mathrm{2}}})$

11:

compute d_LR=LowRankDiag(Z_kΔ_k,Z_k) using Algorithm 2 and d_s=diagQ_s

12:

estimate diagonal of approximate posterior covariance matrix ${\mathit{d}}_{k,\mathit{\lambda }}={\mathit{\lambda }}^{-\mathrm{2}}{\mathit{d}}_{\mathrm{s}}-{\mathit{d}}_{\mathrm{LR}}$

13:

return  MAP estimate s_k,λ and variance estimate d_k,λ

In the following, we provide some details regarding the estimation of the posterior variance in the known mean case. This material has previously appeared in Chung et al. (2018, Sect. 4.1), but we provide a brief description here for completeness. An alternative expression for the posterior covariance is

$${\displaystyle }{\mathbf{Q}}_{\mathrm{post}}={\mathbf{Q}}_{\mathrm{s}}({\mathit{\lambda }}^{\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}+{\mathbf{Q}}_{\mathrm{s}}{\mathbf{H}}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathbf{HQ}}_{\mathrm{s}}{)}^{-\mathrm{1}}{\mathbf{Q}}_{\mathrm{s}},$$

which is obtained by factoring out Q_s. This expression is not computationally feasible for large inverse problems but can be approximated using the outputs of the generalized Golub–Kahan bidiagonalization (described in Sect. 3.1.1). After k steps of the generalized Golub–Kahan bidiagonalization approach, we have matrices ${\mathbf{U}}_{k+\mathrm{1}},{\mathbf{V}}_k$, and B_k. Let ${\mathbf{B}}_k^{\top }{\mathbf{B}}_k={\mathbf{W}}_k{\mathbf{\Theta }}_k{\mathbf{W}}_k^{\top }$ be the eigenvalue decomposition with eigenvalues ${\mathit{\theta }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\theta }}_k$. Next, we compute the matrix Z_k=Q_sV_kW_k and the diagonal matrix

$${\displaystyle }{\mathbf{\Delta }}_k={\mathit{\lambda }}^{-\mathrm{2}}\left[\begin{array}{c}\frac{{\mathit{\theta }}_{\mathrm{1}}}{{\mathit{\theta }}_{\mathrm{1}}+{\mathit{\lambda }}^{\mathrm{2}}}\\ & \mathrm{\ddots }\\ & & & \frac{{\mathit{\theta }}_k}{{\mathit{\theta }}_k+{\mathit{\lambda }}^{\mathrm{2}}}\end{array}\right]\in {\mathbb{R}}^{k\times k}.$$

Then we can approximate the posterior covariance matrix as

$${\displaystyle }{\stackrel{\mathrm{̃}}{\mathbf{Q}}}_{\mathrm{post}}={\mathbf{Q}}_{\mathrm{s}}({\mathit{\lambda }}^{\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}+{\mathbf{Z}}_k{\mathbf{\Theta }}_k{\mathbf{Z}}_k^{\top }{)}^{-\mathrm{1}}{\mathbf{Q}}_{\mathrm{s}}={\mathit{\lambda }}^{-\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}-{\mathbf{Z}}_k{\mathbf{\Delta }}_k{\mathbf{Z}}_k^{\top },$$

providing an efficient representation of the posterior covariance matrix (as a low-rank perturbation of the prior covariance matrix) that can be used for efficient uncertainty quantification. The accuracy of the approximate posterior covariance matrix and of the resulting posterior distribution can be monitored using the information available from the generalized Golub–Kahan bidiagonalization (Saibaba et al., 2020). Note that in the hybrid approach, we estimate the regularization parameters typically at every iteration; the results in Saibaba et al. (2020) are to be applied after the iterations have been terminated and the regularization parameter has been estimated. The uncertainty estimates therefore depend on the value of the regularization parameter. In general, the approximate posterior variance overestimates the uncertainty and decreases monotonically with the number of iterations k.

To visualize the uncertainty, it is common to compute the diagonals of the posterior covariance (known as the posterior variance). The diagonals of Q_s are typically known analytically; the diagonals of the low-rank term ${\mathbf{Z}}_k{\mathbf{\Delta }}_k{\mathbf{Z}}_k^{\top }$ can be computed efficiently using Algorithm 2 with input Y=Z_kΔ_k and Z=Z_k. An important point worth emphasizing is that the approximation to the posterior covariance need not be computed explicitly. More precisely, in addition to storing the information required for storing Q_s, we only need to store nk+k additional entries corresponding to the matrices Z_k and Δ_k.

In addition, one can also use genHyBRs to compute the posterior variance of the sum of the fluxes (or analogously the variance of the mean). To do so, let 1 denote an n×1 vector of ones and multiply the components of ${\stackrel{\mathrm{̃}}{\mathbf{Q}}}_{\mathrm{post}}$ as

$${\displaystyle }{\mathbf{1}}^{\top }{\stackrel{\mathrm{̃}}{\mathbf{Q}}}_{\mathrm{post}}\mathbf{1}={\mathit{\lambda }}^{-\mathrm{2}}\left({\mathbf{1}}^{\top }{\mathbf{Q}}_{\mathrm{s}}\mathbf{1}\right)-\left({\mathbf{1}}^{\top }{\mathbf{Z}}_k\right){\mathbf{\Delta }}_k({\mathbf{1}}^{\top }{\mathbf{Z}}_k{)}^{\top }.$$

Several existing studies (Yadav and Michalak, 2013; Miller et al., 2020) describe how to efficiently compute 1^⊤Q_s1 using Kronecker products.

Note that in previous works (Chung et al., 2018; Saibaba et al., 2020), we found that additional re-orthogonalization of the generalized Golub–Kahan basis vectors yielded more accurate results, so we perform them in the numerical experiments.

Algorithm 2 Compute the diagonals of the low-rank term YZ^⊤. Call as [d]=LowRankDiag(Y,Z).

Require:

Matrices $\mathbf{Y},\mathbf{Z}\in {\mathbb{R}}^{n\times k}$ defining the outer product YZ^⊤

1:

for  $i=\mathrm{1},\mathrm{\dots },n$ do

2:

$d_i={\sum }_{j=\mathrm{1}}^k{\mathbf{Y}}_{ij}{\mathbf{Z}}_{ij}$

3:

end for

4:

return  vector d∈ℝⁿ containing the diagonals of YZ^⊤

3.3 Hierarchical Gaussian priors: reformulation for mean estimation

Next we describe genHyBR methods for AIM with unknown means as described in Sect. 2 with assumptions given in Eq. (5). First, we reformulate the problem for simultaneous estimation of the surface fluxes in s and the covariate parameters in β. Then, we describe how to use genHyBR methods for computing the corresponding MAP estimate Eq. (7) and for subsequent UQ. We refer to this approach as genHyBRmean with UQ, and since the derivations follow those in Sect. 3.1 and 3.2, specific details have been relegated to Appendix B. However, we would like to emphasize that this derivation is new and a contribution of this work.

Note that optimization problem Eq. (7) can be rewritten in standard least squares form,

$${\displaystyle }{\mathit{\gamma }}_{\mathrm{post}}{\displaystyle }=\underset{\mathit{\gamma }=[{\mathit{s}}^{\top },{\mathit{\beta }}^{\top }{]}^{\top }}{\mathrm{argmin}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{∥\left[\begin{array}{cc}{\mathbf{L}}_{\mathbf{R}}\mathbf{H}& \mathbf{0}\\ \mathit{\lambda }{\mathbf{L}}_{\mathit{s}}& -\mathit{\lambda }{\mathbf{L}}_{\mathit{s}}\mathbf{X}\\ \mathbf{0}& {\mathit{\lambda }}_{\mathit{\beta }}{\mathbf{L}}_{\mathit{\beta }}\end{array}\right]\left[\begin{array}{c}\mathit{s}\\ \mathit{\beta }\end{array}\right]-\left[\begin{array}{c}{\mathbf{L}}_{\mathbf{R}}\mathit{z}\\ \mathbf{0}\\ {\mathit{\lambda }}_{\mathit{\beta }}{\mathbf{L}}_{\mathit{\beta }}{\mathit{\mu }}_{\mathit{\beta }}\end{array}\right]∥}_{\mathrm{2}}^{\mathrm{2}},$$

if symmetric decompositions ${\mathbf{R}}^{-\mathrm{1}}={\mathbf{L}}_{\mathbf{R}}^{\top }{\mathbf{L}}_{\mathbf{R}},{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}={\mathbf{L}}_{\mathrm{s}}^{\top }{\mathbf{L}}_{\mathrm{s}},$ and ${\mathbf{Q}}_{\mathit{\beta }}^{-\mathrm{1}}={\mathbf{L}}_{\mathit{\beta }}^{\top }{\mathbf{L}}_{\mathit{\beta }}$ are available. Since computing L_s can be computationally unfeasible (e.g., for spatiotemporal covariance matrices where n is large), we propose a similar change in variables to avoid L_s,

$$\begin{array}{}\text{(17)}& \stackrel{\mathrm{̃}}{\mathit{s}}\leftarrow \mathit{s}-\mathbf{X}{\mathit{\mu }}_{\mathit{\beta }},\stackrel{\mathrm{̃}}{\mathit{\beta }}\leftarrow \mathit{\beta }-{\mathit{\mu }}_{\mathit{\beta }},\stackrel{\mathrm{̃}}{\mathit{z}}\leftarrow \mathit{z}-\mathbf{HX}{\mathit{\mu }}_{\mathit{\beta }}.\end{array}$$

We define the concatenated vector $\stackrel{\mathrm{̃}}{\mathit{\gamma }}=[{\stackrel{\mathrm{̃}}{\mathit{s}}}^{\top },{\stackrel{\mathrm{̃}}{\mathit{\beta }}}^{\top }{]}^{\top }$ and let $\mathbf{K}=\left[\begin{array}{cc}\mathbf{H}& \mathbf{0}\end{array}\right]\in {\mathbb{R}}^{m\times (n+p)}$, where $\mathbf{0}\in {\mathbb{R}}^{m\times p}$ is included to ensure proper matrix multiplication. Then, optimization problem Eq. (7) can be written as

$$\begin{array}{}\text{(18)}& \underset{\stackrel{\mathrm{̃}}{\mathit{\gamma }}\in {\mathbb{R}}^{n+p}}{min}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{∥\mathbf{K}\stackrel{\mathrm{̃}}{\mathit{\gamma }}-\stackrel{\mathrm{̃}}{\mathit{z}}∥}_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}{∥\stackrel{\mathrm{̃}}{\mathit{\gamma }}∥}_{{\mathbf{Q}}^{-\mathrm{1}}}^{\mathrm{2}}\end{array}$$

where

$$\begin{array}{}\text{(19)}& \begin{array}{rl}\mathbf{Q}& =\left[\begin{array}{cc}{\mathbf{Q}}_{\mathrm{s}}+\frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}{\mathbf{XQ}}_{\mathit{\beta }}{\mathbf{X}}^{\top }& \frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}{\mathbf{XQ}}_{\mathit{\beta }}\\ \frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}({\mathbf{XQ}}_{\mathit{\beta }}{)}^{\top }& \frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}{\mathbf{Q}}_{\mathit{\beta }}\end{array}\right]\\ & =\left[\begin{array}{cc}{\mathbf{Q}}_{\mathrm{s}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right]+{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}}\left[\begin{array}{c}\mathbf{X}\\ \mathbf{I}\end{array}\right]{\mathbf{Q}}_{\mathit{\beta }}\left[\begin{array}{cc}{\mathbf{X}}^{\top }& \mathbf{I}\end{array}\right].\end{array}\end{array}$$

Derivations are provided in Appendix B1. Note that, in practice, neither of the matrices H or K need to be formed explicitly since we only need access to matrix–vector products with these matrices and their transposes. Also, we do not explicitly construct Q but instead provide an efficient way to form matrix–vector products with Q.

In summary, to handle AIM with an unknown mean, the genHyBR method can be used to solve Eq. (18) (as described in Sect. 3.1 with Q instead of Q_s, K instead of H, and $\stackrel{\mathrm{̃}}{\mathit{z}}$ instead of b) to efficiently obtain the solution ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{k,\mathit{\lambda }}=[{\stackrel{\mathrm{̃}}{\mathit{s}}}_{k,\mathit{\lambda }}^{\top },{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{k,\mathit{\lambda }}^{\top }{]}^{\top }$. Then, we recover the MAP estimate for s and β as

$$\begin{array}{}\text{(20)}& {\mathit{\gamma }}_{\mathrm{post}}\approx {\mathit{\gamma }}_{k,\mathit{\lambda }}:=\left[\begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{s}}}_{k,\mathit{\lambda }}+\mathbf{X}{\mathit{\mu }}_{\mathit{\beta }}\\ {\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{k,\mathit{\lambda }}+{\mathit{\mu }}_{\mathit{\beta }}\end{array}\right].\end{array}$$

Similarly to the fixed mean case, we can efficiently approximate the posterior covariance matrix and its diagonals using elements of the generalized Golub–Kahan bidiagonalization algorithm, as described in Appendix B3. We note that efficient UQ approaches for the unknown mean case were considered in Saibaba and Kitanidis (2015), but our approach differs in that we reuse information contained in the subspaces generated during the iterative method rather than randomization techniques, making these derivations straightforward and the approaches widely applicable.

4 Numerical results

We evaluate the inverse modeling algorithms described in this paper using the two case studies described in Sect. 4.1. For the numerical experiments, we describe the two methods we test.

• genHyBRs refers to the fixed mean case and involves solving the optimization problem (4) using the approach described in Sect. 3.1.1 and 3.1.2.

• genHyBRmean refers to the unknown mean case and involves solving the optimization problem (18) using the approach described in Sect. 3.3.

For comparison, we use a direct inversion method, which solves Eq. (9) using MATLAB's “backslash” operator. Numerical experiments presented here were obtained using MATLAB on a compute server with four Intel 15-core 2.8 GHz processors and 1 TB of RAM.

4.1 Overview of the case studies

We explore two case studies on estimating CO₂ fluxes across North America using observations from NASA's OCO-2 satellite. In the first case study, we estimate CO₂ fluxes for 6 weeks (late June through July 2015), an inverse problem that is small enough to estimate using the direct method. The second case study using 1 year of observations (September 2014–August 2015) is too large to estimate directly on many or most computer systems. The goal of these experiments is to demonstrate the performance of the generalized hybrid methods for solving the inverse problem with automatic parameter selection.

The case studies explored here parallel those in Miller et al. (2020) and Liu et al. (2021). We provide an overview of these case studies but refer to Miller et al. (2020) for additional detail on the specific setup. Both of the case studies use synthetic OCO-2 observations that are generated using CO₂ fluxes from the NOAA's CarbonTracker product (version 2019b). As a result of this setup, the true CO₂ fluxes (s) are known, making it easier to evaluate the accuracy of the algorithms tested here. All atmospheric transport simulations are from the Weather Research and Forecasting (WRF) Stochastic Time-Inverted Lagrangian Transport Model (STILT) modeling system (Lin et al., 2003; Nehrkorn et al., 2010). These simulations were generated as part of the NOAA's CarbonTracker-Lagrange program (Hu et al., 2019; Miller et al., 2020). Note that the WRF-STILT outputs can be used to explicitly construct H, making it straightforward to calculate the direct solution to the inverse problem in the 6-week case study. By contrast, many inverse modeling studies use the adjoint of an Eulerian model. These modeling frameworks rarely produce an explicit H but instead output the product of H or H^⊤ and a vector. Though we use WRF-STILT for the case studies presented here, the genHyBR algorithms could also be paired with the adjoint of an Eulerian model.

The goal is to estimate CO₂ fluxes at a 3-hourly temporal resolution and a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ latitude–longitude spatial resolution. Using this modeling framework, synthetic observations were obtained as in Eq. (1) by adding white Gaussian noise (representing measurement and model errors) to the output from the atmospheric transport model. Note that WRF-STILT footprints from CarbonTracker-Lagrange are available at 2 s intervals along the OCO-2 flight track, meaning that there are fewer footprints available than there are observations in the original OCO-2 data files. Indeed, many inverse modeling studies to date have thinned or averaged the OCO-2 observations before assimilating them within an inverse model (e.g., Peiro et al., 2022). We generate the synthetic observations by multiplying CO₂ fluxes from CarbonTracker by these footprints, described in greater detail in Miller et al. (2020) and Liu et al. (2021). In total, there are $m=\mathrm{1.92}\times {\mathrm{10}}^{\mathrm{4}}$ synthetic observations and $n=\mathrm{1.06}\times {\mathrm{10}}^{\mathrm{6}}$ unknown CO₂ fluxes in the 6-week case study. By contrast, for the much larger 1-year case study, there are $m=\mathrm{9.9}\times {\mathrm{10}}^{\mathrm{4}}$ synthetic observations and $n=\mathrm{9.4}\times {\mathrm{10}}^{\mathrm{6}}$ unknown CO₂ fluxes to be estimated.

The noise covariance matrix is structured as R=σ²I, where σ² represents the noise variance. In this case, the discrepancy principle formula simplifies to

$${\mathcal{D}}_{\mathrm{full}}\left(\mathit{\lambda }\right)=\Vert \mathbf{H}{\mathit{s}}_{\mathit{\lambda }}-\mathit{z}{\Vert }_{\mathrm{2}}^{\mathrm{2}}\le \mathit{\tau }m{\mathit{\sigma }}^{\mathrm{2}}.$$

Note that the DP approach requires a priori knowledge or an estimate of the noise variance σ². In Miller et al. (2020), σ=2 was used, which leads to a relatively large amount of noise in the observations. We test different values of σ, corresponding to different noise levels (referred to as nlevel), as shown in Table 1. More specifically, let n be a realization of the Gaussian process $\sim \mathcal{N}(\mathbf{0},\mathbf{I})$; then the amount of noise added in Eq. (1) is ϵ=σn, where $\mathit{\sigma }=\mathtt{\text{nlevel}}\cdot \frac{\Vert \mathbf{H}\mathit{s}{\Vert }_{\mathrm{2}}}{\Vert \mathit{n}{\Vert }_{\mathrm{2}}}$. We note that some of the considered noise levels, although very high compared to examples in the inverse problems literature, are lower than typically observed in practice for atmospheric inverse problems. However, recent studies by Miller et al. (2018), O'Dell et al. (2018), Crowell et al. (2019), and Miller and Michalak (2020) show that errors in OCO-2 observations have been gradually decreasing with regular improvements in the satellite retrieval algorithms and bias corrections. Some of the values in Table 1 are low, even considering these recent improvements. With that said, these values are aspirational and may become more realistic in the future as observational and atmospheric modeling errors decline. Furthermore, they provide an opportunity to explore the behavior of the proposed inverse modeling algorithms at many different error levels.

Table 1 Noise level and corresponding noise standard deviation σ used in the 6-week case study experiments.

Figure 2 Sparsity pattern of the prior covariance matrices Q_t and Q_g for the 6-week case study.

Next, we describe the prior used in both case studies. The prior flux estimate is set to a constant value for the case studies explored here. As a result of this setup, the prior flux estimate does not contain any spatiotemporal patterns, and the patterns in the posterior fluxes solely reflect the information content of the atmospheric observations. For the cases with a fixed mean (genHyBRs), we set μ=0, as has been done in several existing studies on inverse modeling algorithms (Rodgers, 2000; Chung and Saibaba, 2017). For the unknown mean case (genHyBRmean), the columns of X contain ones and zeros, denoting fluxes in a given time period. Specifically, for the 6-week case study, X has eight columns, and each column corresponds to a different 3-hourly time period of the day. A given column of X contains values of one for all flux elements that correspond to a given 3-hourly time of day and zero for all other elements. CO₂ fluxes have a large diurnal cycle, and this setup accounts for the fact that CO₂ fluxes at different times of day will have a different mean. In the 1-year case study, X has 12 columns, and each column corresponds to fluxes in a different month of the year, a setup identical to that used in Miller et al. (2020).

For the prior covariance matrix of unknown fluxes, ${\mathbf{Q}}_{\mathrm{s}}={\mathbf{Q}}_{\mathrm{t}}\otimes {\mathbf{Q}}_g$, where Q_t represents the temporal covariance and Q_g represents the spatial covariance in the fluxes. The symbol ⊗ denotes the Kronecker product. We use a spherical covariance model for the spatial and temporal covariance. A spherical model is ideal because it decays to zero at the correlation length or time, and the resulting matrices are usually sparse:

$$\begin{array}{}\text{(21)}& {\displaystyle }{k}_{\mathrm{t}}(d_{\mathrm{t}};{\mathit{\theta }}_{\mathrm{t}})& {\displaystyle }=\left\{\begin{array}{ll}\mathrm{1}-\frac{\mathrm{3}}{\mathrm{2}}\left(\frac{d_{\mathrm{t}}}{{\mathit{\theta }}_{\mathrm{t}}}\right)+\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{d_{\mathrm{t}}}{{\mathit{\theta }}_{\mathrm{t}}}\right)}^{\mathrm{3}}& \text{ if }{d}_{\mathrm{t}}\le {\mathit{\theta }}_{\mathrm{t}},\\ \mathrm{0}& \text{ if }{d}_{\mathrm{t}}>{\mathit{\theta }}_{\mathrm{t}},\end{array}\right.\text{(22)}& {\displaystyle }{k}_g(d_g;{\mathit{\theta }}_g)& {\displaystyle }=\left\{\begin{array}{ll}\mathrm{1}-\frac{\mathrm{3}}{\mathrm{2}}\left(\frac{d_g}{{\mathit{\theta }}_g}\right)+\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{d_g}{{\mathit{\theta }}_g}\right)}^{\mathrm{3}}& \text{ if }{d}_g\le {\mathit{\theta }}_g,\\ \mathrm{0}& \text{ if }{d}_g>{\mathit{\theta }}_g,\end{array}\right.\end{array}$$

where d_t is the temporal difference, d_g is the spherical distance, and θ_t,θ_g represent the decorrelation time, which has units of d, and decorrelation length, which has units of km, respectively. For the 6-week case study, we set θ_t=9.854 d and θ_g=555.42 km, as in Miller et al. (2020). The sparsity patterns of these covariance matrices are provided in Fig. 2. We further set the diagonal elements of Q_s equal to one, and the regularization parameter in the inverse model (λ) will ultimately scale Q_s (e.g., Eq. 2). For the 1-year study, we use slightly different parameters, as listed in full detail in Miller et al. (2020). Notably, the variance is different in each month to better capture the impact of seasonal changes in the variability of CO₂ fluxes. Note that, for the fixed mean case, the covariance matrices are sparse and can be efficiently represented in factored form. However, the hybrid approaches proposed here can handle much more complicated cases (see, e.g., Chung and Saibaba, 2017; Chung et al., 2018).

For the unknown mean case, the covariance matrix Q_β is set to be the identity matrix, and α=10 is used in the numerical experiments. We experimented with various choices for α and observed consistently good results with α=10. In all of the numerical experiments, we use the DP approach to select the regularization parameter within genHyBR methods with τ=1.

In subsequent discussion of the case study results, we provide relative reconstruction error norms computed as $\left|\right|{\mathit{s}}_{k,\mathit{\lambda }}-\mathit{s}|{|}_{\mathrm{2}}/|\left|\mathit{s}\right|{|}_{\mathrm{2}}$, where s denotes the true fluxes and s_k,λ contains the reconstructed spatiotemporal fluxes at the kth iteration.

4.2 Results of the case studies

Both the genHyBRs and genHyBRmean methods converge quickly and yield accurate estimates of the CO₂ fluxes relative to other inverse modeling methods. For the 6-week case study, Fig. 3 shows the relative reconstruction error norms for both genHyBRs and genHyBRmean for three different noise levels and for different options for selecting the regularization parameter. DP corresponds to the discrepancy principle and is an automatic approach that depends on the data and noise level. The optimal regularization parameter, which corresponds to selecting the regularization parameter at each iteration that minimizes the reconstruction error, is provided for comparison, although it is not obtainable in practice. All of the plots with “none” correspond to λ=0 and show semiconvergent behavior, which is revealed in the “U” shape of the relative error plot. That is, the relative reconstruction error norms decrease in early iterations but increase in later iterations due to noise contamination in the reconstructions.

Figure 3 Six-week case study: relative reconstruction error norms per iteration of genHyBRs and genHyBRmean with 5 %, 10 %, and 50 % noise levels. We compare results for the optimal regularization parameter, the automatically selected DP parameter, λ=0, and the L-BFGS with a fixed regularization parameter (λ=1). The filled circle indicates the stopping iteration for the genHyBR methods with DP. In addition, the horizontal line denotes to the relative error for the reconstruction obtained using the direct method.

Table 2 Six-week case study: for various noise levels, we provide comparisons of genHyBRs and genHyBRmean (with a DP-selected regularization parameter) to standard direct and iterative methods (with a fixed regularization parameter). We provide the number of iterations, the CPU timing in seconds, and the relative reconstruction error norms for the computed spatiotemporal fluxes denoted by Δs. Note that, when computing the reconstructions and approximating the posterior variance, additional reorthogonalization is performed, which explains the slight difference in the number of iterations and run times.

We observe that, for all noise levels, genHyBR methods with regularization parameter estimation (genHyBRs--opt and genHyBRs--dp) result in reconstruction error norms that decrease and flatten, thereby overcoming the semiconvergent behavior of genHyBR methods with no regularization (genHyBRs--none). For reference, we mark with a horizontal line the relative reconstruction error norm for a direct reconstruction. Recall that the direct method and the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method typically require λ to be fixed in advance, and a poor choice of λ can yield poor reconstructions of the CO₂ fluxes. Using the optimal regularization parameter computed from genHyBRmean (these values are provided in Table 2), we show that a good reconstruction can be obtained with the direct method if a good regularization parameter is available, but obtaining this result may require careful and expensive tuning.

Figure 4 Six-week case study: reconstructed fluxes, averaged over 6 weeks, are provided for genHyBRs and genHyBRmean for the automatically selected DP parameter and for λ=0. The true average fluxes and the reconstruction using a direct method with the optimal regularization parameter computed from genHyBRmean are provided for comparison. These results correspond to the 50 % noise level, and relative error norms of average fluxes over 6 weeks are provided in the titles.

Figure 5 Estimated uncertainties for the total CO₂ flux from the continental US over the 6-week case study (50 % noise levels). The figure compares posterior uncertainties calculated using the direct approach, using GenHyBR and using the low-rank approach described in Miller et al. (2020). Solid lines show uncertainties estimated for simulations with a fixed mean and dashed lines simulations with an unknown mean. In addition, we use the same regularization parameter (λ) for the estimates shown here to make the convergence behavior of the different approaches easier to compare (λ=0.13 for the fixed mean and λ=0.15 for the unknown mean, as in Table 2). GenHyBR reaches uncertainty estimates that are close to the direct estimate more quickly than the low-rank approach, particularly for the unknown mean case. In addition, GenHyBR uses forward and adjoint model runs previously generated as part of the best estimate calculations, yielding large computational savings over the low-rank approach, which requires new model runs.

We also find that the algorithm that simultaneously estimates the mean (genHyBRmean) yields lower errors than the algorithm with a fixed mean (genHyBRs). Notably, the difference in performance between these two algorithms grows as the noise level increases. In other words, the comparative advantage of genHyBRmean is even larger at higher noise levels. This result implies that mean estimation becomes critical for problems with large noise levels.

For comparison, we show the convergence behavior of the L-BFGS algorithm that is commonly used in AIM (Fig. 3). One line (with diamond symbols) shows results using L-BFGS to minimize the inverse model with an unknown mean (with λ=1), and another line (shown with plus symbols) shows a variant of L-BFGS with a data transformation to speed convergence (also with λ=1). The specific approach used here is described in detail in Miller et al. (2020). Both algorithms converge more slowly than the genHyBR algorithms. This result is significant because the main computational cost per iteration (i.e., one matrix–vector multiplication by H and its adjoint) is the same among the algorithms in Fig. 3. In addition, the relative errors increase for L-BFGS at later iterations, though the inverse modeling optimization function continues to decrease at these iterations; the inverse model minimizes the errors with respect to the observations (z) but is not guaranteed to minimize errors with respect to the true fluxes (s).

The estimated fluxes using genHyBR also exhibit spatial patterns that mirror fluxes estimated using a direct (e.g., analytical) approach. Maps of CO₂ fluxes, averaged over 6 weeks and corresponding to a 50 % noise level, are shown in Fig. 4. We provide the true average flux, the genHyBRs and genHyBRmean reconstructions for various parameter choices, and the direct method reconstruction using the optimal regularization parameter from genHyBRmean. Note that these relative reconstruction error values are different than those provided in Fig. 3 because they represent error norms computed on the average image rather than on the native 3 h resolution of the estimated fluxes. Also, to better highlight broad spatial patterns, the color map has been constrained so that average flux estimates above 2 are set to 2 and estimates below −5 are set to −5.

Perhaps surprisingly, the genHyBR algorithms require less computing time than other algorithms tested, including the direct or analytical method. Table 2 displays the measured turnaround time for each case, along with the number of iterations and the relative reconstruction error norms for each method. Compared to the direct method with fixed λ, hybrid methods require less time to compute the estimated CO₂ fluxes. Moreover, since the regularization parameter can be selected automatically, genHyBR methods can obtain results with smaller reconstruction errors.

Figure 6 One-year case study: for 50 % noise level, we provide relative reconstruction error norms per iteration of genHyBRs and genHyBRmean and compare results for the optimal regularization parameter, the automatically selected DP parameter, and λ=0.

Figure 7 One-year case study: reconstructed fluxes, averaged over 1 year, are provided for genHyBRs and genHyBRmean for various parameter choices. These results correspond to the 50 % noise level, and relative error norms of average fluxes over 1 year are provided in the titles.

Finally, we demonstrate the ability to perform UQ for AIM. In the last columns of Table 2, we provide the times needed to compute uncertainties. We remark that the additional time and the difference in the number of iterations can be attributed to the need to perform reorthogonalization of the Krylov basis vectors, which is not as critical for obtaining the MAP estimate. Furthermore, uncertainty estimation for genHyBR does not require any additional forward or adjoint model runs beyond what is required to estimate the fluxes (s). Hence, the additional computing time is modest, given the size of the problem and the ability to obtain solution variance estimates. The estimated posterior variance is also similar across both methods (genHyBRs and genHyBRmean).

GenHyBR also requires relatively few iterations (i.e., Krylov basis vectors) to converge on a reasonable uncertainty estimate. Figure 5 shows the posterior uncertainty in the total continental US CO₂ flux by iteration for each algorithm. This figure compares the direct approach to uncertainty estimation, GenHyBR, and the low-rank approach described in Saibaba and Kitanidis (2015) and Miller et al. (2020). In this figure, we use the same regularization parameter (λ) across all algorithms to make the results from different algorithms more easily comparable to one another. The low-rank approach described in Saibaba and Kitanidis (2015) and Miller et al. (2020) starts at very high values and converges more slowly than the uncertainties estimated using GenHyBR. Furthermore, the low-rank approach requires new forward and adjoint runs for the uncertainty calculations, while GenHyBR uses forward and adjoint runs already generated when calculating the best estimate of the fluxes.

We conclude this section with results for the 1-year case study. Since this case study has approximately 9 times the number of unknown CO₂ fluxes and 5 times the number of observations compared to the 6-week case study, iterative methods are computationally more appealing than direct methods for obtaining reconstructions. The previous case study already explored the behavior of the algorithms at different noise levels, so here we only consider the 50 % noise level, which corresponds to σ=0.4076. Figure 6 provides relative reconstruction error norms for genHyBRs and genHyBRmean. With the regularization parameter automatically selected using DP, both methods result in reconstructions with relative errors smaller than 0.85. Since it is difficult to show spatiotemporal flux reconstructions over the entire year, we provide the annual average of the reconstructed CO₂ fluxes in Fig. 7. Compared to the reconstructions in Fig. 4, these average maps are much smoother. This inability to resolve fine details can be attributed to the significantly fewer observations compared to the number of unknowns in this case study. Nevertheless, these results show that the algorithms described in this study can be scaled to very large inverse problems – problems where the direct method is either computationally prohibitive or time-consuming.

5 Conclusions

This article describes a mathematically advanced iterative method for AIM with large data sets. Specifically, we discuss generalized hybrid methods for inverse models with a fixed prior mean (e.g., Bayesian synthesis inverse modeling) and an unknown prior mean (e.g., geostatistical inverse modeling). We also describe a means of obtaining posterior variance estimates at very little additional computational cost. Compared to standard inverse modeling procedures (e.g., direct and iterative methods), genHyBR methods are computationally cheaper and exhibit faster convergence. One of the main advantages of genHyBR methods is the ability to efficiently and adaptively estimate the regularization parameter during the inversion process, and we described various regularization parameter estimation methods for Tikhonov regularization. Numerical experiments for case studies for 6 weeks and 1 year demonstrate that genHyBR methods provide an efficient, flexible, robust, and automatic approach for AIM with very large spatiotemporal fluxes. Furthermore, since these methods only require forward and adjoint model evaluations, these methods can be paired with different types of atmospheric transport models.

Appendix A: Regularization parameter estimation methods for genHyBR

One of the main advantages of hybrid projection methods is the ability to adaptively and automatically estimate the regularization parameter during the iterative process. We described the discrepancy principle (DP), but other common parameter estimation techniques in the context of hybrid projection methods include the generalized cross-validation (GCV) method and the weighted GCV (WGCV) method (Chung et al., 2008; Renaut et al., 2017). A summary of methods with respective functions used to compute the parameters based on the original problem and the projected problem are summarized in Table A1, where for simplicity we have assumed that R=σ²I. We have used the notation ${\mathbf{B}}_{k,\mathit{\lambda }}^{†}=({\mathbf{B}}_k^{\top }{\mathbf{B}}_k+{\mathit{\lambda }}^{\mathrm{2}}\mathbf{I}{)}^{-\mathrm{1}}{\mathbf{B}}_k^{\top }$ for given λ>0, and y_k,λ is the solution to the projected problem (14).

Table A1 Regularization parameter selection methods for use within genHyBR methods.

More specifically, DP selects the largest parameter value λ for which 𝒟_full(λ)≤τmσ², where τ≥1 is a user-defined parameter. Note that mσ² is the expected value of $\Vert \mathit{ϵ}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}$. For the projected problem, we choose the largest λ such that ${\mathcal{D}}_{\mathrm{proj}}\left(\mathit{\lambda }\right)\le \mathit{\tau }m{\mathit{\sigma }}^{\mathrm{2}}.$ The WGCV method selects λ by minimizing the objective function 𝒢_full(λ;ω). Note that, if ω=1, then WGCV becomes GCV. In the projected problem, we minimize 𝒢_proj(λ;ω) at each iteration. The parameter ω can be chosen automatically, as described in Chung et al. (2008) and Renaut et al. (2017). Note that the DP approach requires a priori knowledge of the noise variance σ², whereas GCV approaches do not require prior knowledge about the noise level.

Appendix B: Extension to the unknown mean: hierarchical Bayes

In this section, we provide details for the derivation of genHyBR methods for AIM with an unknown mean. We begin in Appendix B1 with the problem reformulation to simultaneously estimate the unknown fluxes and the covariate parameters. Then in Sect. B2 we provide the details of the genHyBR approach for the unknown mean case, which closely follows the derivation in Sect. 3.1.1. Finally, in Appendix B3 we show how to approximate the posterior variance using the generalized Golub–Kahan bidiagonalization.

B1 Reformulation for simultaneous estimation

In order to apply genHyBR methods to the unknown mean estimation problem, we first reformulate the MAP estimate from Eqs. (7) to (18) as follows. For the data fit term, consider

$$\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert \mathbf{H}\mathit{s}-\mathit{z}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{∥\left[\begin{array}{cc}\mathbf{H}& \mathbf{0}\end{array}\right]\left[\begin{array}{c}\mathit{s}\\ \mathit{\beta }\end{array}\right]-\mathit{z}∥}_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{∥\mathbf{K}\left[\begin{array}{c}\stackrel{\mathrm{̃}}{\mathit{s}}+\mathbf{X}{\mathit{\mu }}_{\mathit{\beta }}\\ \stackrel{\mathrm{̃}}{\mathit{\beta }}+{\mathit{\mu }}_{\mathit{\beta }}\end{array}\right]-\mathit{z}∥}_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{∥\mathbf{K}\left[\begin{array}{c}\stackrel{\mathrm{̃}}{\mathit{s}}\\ \stackrel{\mathrm{̃}}{\mathit{\beta }}\end{array}\right]+\mathbf{HX}{\mathit{\mu }}_{\mathit{\beta }}-\mathit{z}∥}_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{∥\mathbf{K}\left[\begin{array}{c}\stackrel{\mathrm{̃}}{\mathit{s}}\\ \stackrel{\mathrm{̃}}{\mathit{\beta }}\end{array}\right]-\stackrel{\mathrm{̃}}{\mathit{z}}∥}_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}},\end{array}$$

and for the regularization terms in Eq. (7), we have

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \mathit{s}-\mathbf{X}\mathit{\beta }{\Vert }_{{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\lambda }}_{\mathit{\beta }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \mathit{\beta }-{\mathit{\mu }}_{\mathit{\beta }}{\Vert }_{{\mathbf{Q}}_{\mathit{\beta }}^{-\mathrm{1}}}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \stackrel{\mathrm{̃}}{\mathit{s}}-\mathbf{X}\stackrel{\mathrm{̃}}{\mathit{\beta }}{\Vert }_{{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\lambda }}_{\mathit{\beta }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \stackrel{\mathrm{̃}}{\mathit{\beta }}{\Vert }_{{\mathbf{Q}}_{\mathit{\beta }}^{-\mathrm{1}}}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\lambda }}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{s}}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\stackrel{\mathrm{̃}}{\mathit{s}}-\mathrm{2}{\mathit{\lambda }}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{s}}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}\stackrel{\mathrm{̃}}{\mathit{\beta }}\\ {\displaystyle }& {\displaystyle }+{\stackrel{\mathrm{̃}}{\mathit{\beta }}}^{\top }({\mathit{\lambda }}^{\mathrm{2}}{\mathbf{X}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}+{\mathit{\lambda }}_{\mathit{\beta }}^{\mathrm{2}}{\mathbf{Q}}_{\mathit{\beta }}^{-\mathrm{1}})\stackrel{\mathrm{̃}}{\mathit{\beta }})\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\lambda }}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{s}}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\stackrel{\mathrm{̃}}{\mathit{s}}-\mathrm{2}{\mathit{\lambda }}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{s}}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}\stackrel{\mathrm{̃}}{\mathit{\beta }}\\ {\displaystyle }& {\displaystyle }+{\stackrel{\mathrm{̃}}{\mathit{\beta }}}^{\top }({\mathit{\lambda }}^{\mathrm{2}}{\mathbf{X}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}+(\mathit{\alpha }\mathit{\lambda }{)}^{\mathrm{2}}{\mathbf{Q}}_{\mathit{\beta }}^{-\mathrm{1}}\left)\stackrel{\mathrm{̃}}{\mathit{\beta }}\right)\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\begin{array}{cc}{\stackrel{\mathrm{̃}}{\mathit{s}}}^{\top }& {\stackrel{\mathrm{̃}}{\mathit{\beta }}}^{\top }\end{array}\right]\\ {\displaystyle }& {\displaystyle }\left[\begin{array}{cc}{\mathit{\lambda }}^{\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}& -{\mathit{\lambda }}^{\mathrm{2}}{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}\\ -{\mathit{\lambda }}^{\mathrm{2}}{\mathbf{X}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}& {\mathit{\lambda }}^{\mathrm{2}}{\mathbf{X}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}+(\mathit{\alpha }\mathit{\lambda }{)}^{\mathrm{2}}{\mathbf{Q}}_{\mathit{\beta }}^{-\mathrm{1}}\end{array}\right]\left[\begin{array}{c}\stackrel{\mathrm{̃}}{\mathit{s}}\\ \stackrel{\mathrm{̃}}{\mathit{\beta }}\end{array}\right]\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}\left[\begin{array}{cc}{\stackrel{\mathrm{̃}}{\mathit{s}}}^{\top }& {\stackrel{\mathrm{̃}}{\mathit{\beta }}}^{\top }\end{array}\right]\\ {\displaystyle }& {\displaystyle }\underset{={\mathbf{Q}}^{-\mathrm{1}}}{\underbrace{\left[\begin{array}{cc}{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}& -{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}\\ -{\mathbf{X}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}& {\mathbf{X}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}+{\mathit{\alpha }}^{\mathrm{2}}{\mathbf{Q}}_{\mathit{\beta }}^{-\mathrm{1}}\end{array}\right]}}\left[\begin{array}{c}\stackrel{\mathrm{̃}}{\mathit{s}}\\ \stackrel{\mathrm{̃}}{\mathit{\beta }}\end{array}\right]\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{2}}}{\mathrm{2}}}\Vert \stackrel{\mathrm{̃}}{\mathit{\gamma }}{\Vert }_{{\mathbf{Q}}^{-\mathrm{1}}}^{\mathrm{2}}.\end{array}$$

We identify $\stackrel{\mathrm{̃}}{\mathit{\gamma }}=[{\stackrel{\mathrm{̃}}{\mathit{s}}}^{\top },{\stackrel{\mathrm{̃}}{\mathit{\beta }}}^{\top }{]}^{\top }$, and this completes the derivation of Eq. (18).

Next we provide the derivation of the augmented prior covariance matrix Eq. (19). Since we need Q and not Q⁻¹ in genHyBR, we use the formula for the inverse of a 2×2 block matrix:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\left[\begin{array}{c}\mathbf{AB}\\ \mathbf{C}& \mathbf{D}\end{array}\right]}^{-\mathrm{1}}=\\ {\displaystyle }& {\displaystyle }\left[\begin{array}{cc}(\mathbf{A}-{\mathbf{BD}}^{-\mathrm{1}}\mathbf{C}{)}^{-\mathrm{1}}& -(\mathbf{A}-{\mathbf{BD}}^{-\mathbf{1}}\mathbf{C}{)}^{-\mathbf{1}}{\mathbf{BD}}^{-\mathbf{1}}\\ -{\mathbf{D}}^{-\mathrm{1}}\mathbf{C}(\mathbf{A}-{\mathbf{BD}}^{-\mathrm{1}}\mathbf{C}{)}^{-\mathrm{1}}& {\mathbf{D}}^{-\mathrm{1}}\mathbf{C}(\mathbf{A}-{\mathbf{BD}}^{-\mathrm{1}}\mathbf{C}{)}^{-\mathrm{1}}{\mathbf{BD}}^{-\mathrm{1}}+{\mathbf{D}}^{-\mathrm{1}},\end{array}\right]\end{array}$$

which is defined if D and $\mathbf{A}-{\mathbf{BD}}^{-\mathrm{1}}\mathbf{C}$ are invertible. Using the above formula, we get

$$\begin{array}{ll}{\displaystyle }\mathbf{Q}& {\displaystyle }={\left[\begin{array}{cc}{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}& -{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}\\ -{\mathbf{X}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}& {\mathbf{X}}^{\top }{\mathbf{Q}}_{\mathrm{s}}^{-\mathrm{1}}\mathbf{X}+{\mathit{\alpha }}^{\mathrm{2}}{\mathbf{Q}}_{\mathit{\beta }}^{-\mathrm{1}}\end{array}\right]}^{-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }=\left[\begin{array}{cc}{\mathbf{Q}}_{\mathrm{s}}+\frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}{\mathbf{XQ}}_{\mathit{\beta }}{\mathbf{X}}^{\top }& \frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}{\mathbf{XQ}}_{\mathit{\beta }}\\ \frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}{\mathbf{Q}}_{\mathit{\beta }}{\mathbf{X}}^{\top }& \frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}{\mathbf{Q}}_{\mathit{\beta }}\end{array}\right].\end{array}$$

This simplifies to

$$\mathbf{Q}=\left[\begin{array}{cc}{\mathbf{Q}}_{\mathrm{s}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right]+{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}}\left[\begin{array}{c}\mathbf{X}\\ \mathbf{I}\end{array}\right]{\mathbf{Q}}_{\mathit{\beta }}\left[\begin{array}{cc}{\mathbf{X}}^{\top }& \mathbf{I}\end{array}\right],$$

which completes the derivation of Eq. (19).

B2 genHyBR approach for AIM with an unknown mean

In this subsection, we derive the genHyBR approach for the unknown mean case. We initialize ${\mathit{\delta }}_{\mathrm{1}}^K=\Vert \stackrel{\mathrm{̃}}{\mathit{z}}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}$, ${\mathit{u}}_{\mathrm{1}}^K=\stackrel{\mathrm{̃}}{\mathit{z}}/{\mathit{\delta }}_{\mathrm{1}}^K$ and ${\mathit{\gamma }}_{\mathrm{1}}^K{\mathit{v}}_{\mathrm{1}}^K={\mathbf{K}}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathit{u}}_{\mathrm{1}}^K$; then, the kth iteration of the generalized Golub–Kahan bidiagonalization procedure generates vectors ${\mathit{u}}_{k+\mathrm{1}}^K$ and ${\mathit{v}}_{k+\mathrm{1}}^K$ such that

$$\begin{array}{ll}{\displaystyle }{\mathit{\gamma }}_{k+\mathrm{1}}^K{\mathit{u}}_{k+\mathrm{1}}^K=& {\displaystyle }\mathbf{KQ}{\mathit{v}}_k^K-{\mathit{\gamma }}_k^K{\mathit{u}}_k^K,\\ {\displaystyle }{\mathit{\delta }}_{k+\mathrm{1}}^K{\mathit{v}}_{k+\mathrm{1}}^K=& {\displaystyle }{\mathbf{K}}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathit{u}}_{k+\mathrm{1}}^K-{\mathit{\delta }}_{k+\mathrm{1}}^K{\mathit{v}}_k^K,\end{array}$$

where scalars ${\mathit{\gamma }}_k^K,{\mathit{\delta }}_k^K\ge \mathrm{0}$ are computed such that $\Vert {\mathit{u}}_k^K{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}=\Vert {\mathit{v}}_k^K{\Vert }_{\mathbf{Q}}=\mathrm{1}$. At the end of k iterations, we have

$$\begin{array}{ll}{\displaystyle }{\mathbf{B}}_k^K& {\displaystyle }\equiv \left[\begin{array}{c}{\mathit{\gamma }}_{\mathrm{1}}^K\\ {\mathit{\delta }}_{\mathrm{2}}^K& {\mathit{\gamma }}_{\mathrm{2}}^K\\ & {\mathit{\delta }}_{\mathrm{3}}^K& \mathrm{\ddots }\\ & & \mathrm{\ddots }& {\mathit{\gamma }}_k^K\\ & & & {\mathit{\delta }}_{k+\mathrm{1}}^K\end{array}\right],\\ {\displaystyle }& {\displaystyle }{\mathbf{U}}_{k+\mathrm{1}}^K\equiv \left[\begin{array}{c}{\mathit{u}}_{\mathrm{1}}^K,\mathrm{\dots },{\mathit{u}}_{k+\mathrm{1}}^K\end{array}\right],\text{ and }{\mathbf{V}}_k^K\equiv \left[\begin{array}{c}{\mathit{v}}_{\mathrm{1}}^K,\mathrm{\dots },{\mathit{v}}_k^K\end{array}\right].\end{array}$$

The above matrices satisfy the following relations:

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{\mathbf{U}}_{k+\mathrm{1}}^K{\mathit{\delta }}_{\mathrm{1}}^K{\mathit{e}}_{\mathrm{1}}& =\stackrel{\mathrm{̃}}{\mathit{z}},\\ {\mathbf{KQV}}_k^K& ={\mathbf{U}}_{k+\mathrm{1}}^K{\mathbf{B}}_k^K,\\ {\mathbf{K}}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathbf{U}}_{k+\mathrm{1}}^K& ={\mathbf{V}}_k^K({\mathbf{B}}_k^K{)}^{\top }+{\mathit{\gamma }}_{k+\mathrm{1}}^K{\mathit{v}}_{k+\mathrm{1}}^K{\mathit{e}}_{k+\mathrm{1}}^{\top }.\end{array}\end{array}$$

Also, matrices ${\mathbf{U}}_{k+\mathrm{1}}^K$ and ${\mathbf{V}}_k^K$ satisfy the following orthogonality conditions:

$$\begin{array}{}\text{(B2)}& ({\mathbf{U}}_{k+\mathrm{1}}^K{)}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathbf{U}}_{k+\mathrm{1}}^K={\mathbf{I}}_{k+\mathrm{1}}\text{ and }({\mathbf{V}}_k^K{)}^{\top }{\mathbf{QV}}_k^K={\mathbf{I}}_k.\end{array}$$

The columns of the matrix ${\mathbf{V}}_k^K$ form a basis for the Krylov subspace ${\mathcal{K}}_k({\mathbf{K}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathbf{KQ},{\mathbf{K}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\stackrel{\mathrm{̃}}{\mathit{z}})$, which we use to search for approximate solutions.

To obtain the approximate solution, we solve the least squares problem

$$\begin{array}{}\text{(B3)}& \underset{\mathit{y}\in {\mathbb{R}}^k}{min}\Vert {\mathbf{B}}_k^K\mathit{y}-{\mathit{\delta }}_{\mathrm{1}}^K{\mathit{e}}_{\mathrm{1}}{\Vert }_{\mathrm{2}}^{\mathrm{2}}+{\mathit{\lambda }}^{\mathrm{2}}\Vert \mathit{y}{\Vert }_{\mathrm{2}}^{\mathrm{2}}\end{array}$$

to obtain the optimizer ${\mathit{y}}_{k,\mathit{\lambda }}^K$ and to compute the approximate solution ${\mathit{\gamma }}_{k,\mathit{\lambda }}={\mathbf{QV}}_k^K{\mathit{y}}_{k,\mathit{\lambda }}^K$. We can extract the approximations ${\stackrel{\mathrm{̃}}{\mathit{s}}}_{k,\mathit{\lambda }}$ and ${\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{k,\mathit{\lambda }}$ as ${\mathit{\gamma }}_{k,\mathit{\lambda }}=\left[\begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{s}}}_{k,\mathit{\lambda }}\\ {\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{k,\mathit{\lambda }}\end{array}\right]$. To estimate the regularization parameter, we can adapt the techniques in Sect. 3.1.2; for example, using the discrepancy principle, we pick a regularization parameter λ such that

$${\displaystyle }{\mathcal{D}}_{\mathrm{proj}}^K\left(\mathit{\lambda }\right)=\Vert {\mathbf{B}}_k^K{\mathit{y}}_{k,\mathit{\lambda }}^K-{\mathit{\delta }}_{\mathrm{1}}^K{\mathit{e}}_{\mathrm{1}}{\Vert }_{\mathrm{2}}^{\mathrm{2}}\le \mathit{\tau }m,$$

where τ≥1 is a user-defined parameter and m is the expected value of $\Vert \mathit{ϵ}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}^{\mathrm{2}}$. Other parameter selection techniques such as GCV and WGCV can also be adapted to the unknown mean case with similar expressions to Table A1, but we omit the details.

B3 Approximation to the posterior variance

We propose an approximation to the posterior covariance matrix Γ_post corresponding to the posterior distribution $\mathit{\pi }(\mathit{s},\mathit{\beta }|\mathit{z})$ in Eq. (6). First note that, from Eqs. (19) and (8), we get the following expression of the posterior covariance matrix:

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}{\mathbf{\Gamma }}_{\mathrm{post}}& ={\left({\mathit{\lambda }}^{\mathrm{2}}{\mathbf{Q}}^{-\mathrm{1}}+{\mathbf{K}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathbf{K}\right)}^{-\mathrm{1}}\\ & =\mathbf{Q}{\left({\mathit{\lambda }}^{\mathrm{2}}\mathbf{Q}+\mathbf{Q}\left({\mathbf{K}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathbf{K}\right)\mathbf{Q}\right)}^{-\mathrm{1}}\mathbf{Q},\end{array}\end{array}$$

where the last expression is obtained by factoring out Q. Assume that k iterations of the generalized Golub–Kahan bidiagonalization have been performed to solve Eq. (18). Let $({\mathbf{B}}_k^K{)}^{\top }{\mathbf{B}}_k^K={\mathbf{W}}_k^K{\mathbf{\Theta }}_k^K({\mathbf{W}}_k^K{)}^{\top }$ be an eigenvalue decomposition with eigenvalues ${\mathit{\theta }}_{\mathrm{1}}^K,\mathrm{\dots },{\mathit{\theta }}_k^K$ and let ${\mathbf{Z}}_k^K={\mathbf{QV}}_k^K{\mathbf{W}}_k^K$. Then consider the following low-rank approximation:

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}\mathbf{Q}\left({\mathbf{K}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathbf{K}\right)\mathbf{Q}& \approx \mathbf{Q}\left({\mathbf{V}}_k^K\right({\mathbf{B}}_k^K{)}^{\top }{\mathbf{B}}_k^K{\mathbf{V}}_k^{\top })\mathbf{Q}\\ & ={\mathbf{Z}}_k^K{\mathbf{\Theta }}_k^K({\mathbf{Z}}_k^K{)}^{\top }.\end{array}\end{array}$$

Using Eq. (B5) in Eq. (B4), we can define an approximate covariance matrix as

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathbf{\Gamma }}}_{\mathrm{post}}& =\mathbf{Q}({\mathit{\lambda }}^{\mathrm{2}}\mathbf{Q}+{\mathbf{Z}}_k^K{\mathbf{\Theta }}_k^K({\mathbf{Z}}_k^K{)}^{\top }{)}^{-\mathrm{1}}\mathbf{Q}\\ & ={\mathit{\lambda }}^{-\mathrm{2}}\mathbf{Q}-{\mathbf{Z}}_k^K{\mathbf{\Delta }}_k^K({\mathbf{Z}}_k^K{)}^{\top },\end{array}\end{array}$$

where ${\mathbf{\Delta }}_k^K$ is a diagonal matrix,

$${\mathbf{\Delta }}_k^K\equiv {\mathit{\lambda }}^{-\mathrm{2}}\left[\begin{array}{c}\frac{{\mathit{\theta }}_{\mathrm{1}}^K}{{\mathit{\theta }}_{\mathrm{1}}^K+{\mathit{\lambda }}^{\mathrm{2}}}\\ & \mathrm{\ddots }\\ & & \frac{{\mathit{\theta }}_k^K}{{\mathit{\theta }}_k^K+{\mathit{\lambda }}^{\mathrm{2}}}\end{array}\right]\in {\mathbb{R}}^{k\times k}.$$

The last expression was obtained using the Woodbury formula.

Therefore, we have an efficient representation of the approximate posterior covariance matrix as a low-rank perturbation of the prior covariance matrix Q. It is important to note that, as with the prior covariance matrix Q, we do not need to store ${\stackrel{\mathrm{̃}}{\mathbf{\Gamma }}}_{\mathrm{post}}$ explicitly. More precisely, in addition to storing the information required for storing Q, we only need to store nk+k additional entries corresponding to the matrices ${\mathbf{Z}}_k^K$ and ${\mathbf{\Delta }}_k^K$. Furthermore, the error in the low-rank approximation can be analyzed using similar techniques to Saibaba et al. (2020). Similar to the approach described in Sect. 3.2, the posterior variance, which corresponds to the diagonal entries of Γ_post, can be approximated using the diagonal entries of Eq. (B6). First, note that the diagonals of Q are obtained from the diagonals of the block matrices ${\mathbf{Q}}_{\mathrm{s}}+{\mathit{\alpha }}^{-\mathrm{2}}{\mathbf{XQ}}_{\mathit{\beta }}{\mathbf{X}}^{\top }$ and α⁻²Q_β. The diagonals of Q_s and Q_β are typically known analytically. The diagonals of XQ_sX^⊤ and ${\mathbf{Z}}_k^K{\mathbf{\Delta }}_k^K({\mathbf{Z}}_k^K{)}^{\top }$ are easy to compute in 𝒪((n+p)k²) flops since they are low-rank matrices. Therefore, given the approximate representation of the covariance matrix (B6), we can estimate the posterior variance (i.e., the diagonals of the posterior covariance). A complete description of the method is given in Algorithm B1.

Algorithm B1 AIM with unknown mean – genHyBRmean with UQ.

Require:

Matrices K, R and Q, and vector $\stackrel{\mathrm{̃}}{\mathit{z}}$.

1:

{$/*$ Compute MAP estimate $*/$}

2:

initialize ${\mathit{u}}_{\mathrm{1}}^K=\stackrel{\mathrm{̃}}{\mathit{z}}/\Vert \stackrel{\mathrm{̃}}{\mathit{z}}{\Vert }_{{\mathbf{R}}^{-\mathrm{1}}}$

3:

for $j=\mathrm{1},\mathrm{\dots },k$ do

4:

one iteration of generalized Golub–Kahan bidiagonalization to obtain ${\mathbf{B}}_j^K,{\mathbf{U}}_{j+\mathrm{1}}^K,$ and ${\mathbf{V}}_j^K$

5:

estimate regularization parameter λ and compute ${\mathit{y}}_{j,\mathit{\lambda }}^K$ by solving Eq. (B3)

6:

end for

7:

compute $\begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{s}}}_{k,\mathit{\lambda }}\\ {\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{k,\mathit{\lambda }}\end{array}={\mathbf{QV}}_k^K{\mathit{y}}_{k,\mathit{\lambda }}^K$ and ${\mathit{\gamma }}_{k,\mathit{\lambda }}=\left[\begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{s}}}_{k,\mathit{\lambda }}+\mathbf{X}{\mathit{\mu }}_{\mathit{\beta }}\\ {\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{k,\mathit{\lambda }}+{\mathit{\mu }}_{\mathit{\beta }}\end{array}\right]$

8:

{$/*$ Compute the approximation to the posterior variance $*/$}

9:

compute the eigendecomposition

$({\mathbf{B}}_k^K{)}^{\top }{\mathbf{B}}_k^K={\mathbf{W}}_k^K{\mathbf{\Theta }}_k^K({\mathbf{W}}_k^K{)}^{\top }$

10:

compute ${\mathbf{Z}}_k^K={\mathbf{QV}}_k^K{\mathbf{W}}_k^K$ and

${\mathbf{\Delta }}_k^K=\text{diag}(\frac{{\mathit{\theta }}_{\mathrm{1}}^K}{{\mathit{\theta }}_{\mathrm{1}}^K+{\mathit{\lambda }}^{\mathrm{2}}},\mathrm{\dots },\frac{{\mathit{\theta }}_k^K}{{\mathit{\theta }}_k^K+{\mathit{\lambda }}^{\mathrm{2}}})$

11:

compute d_LR=LowRankDiag$({\mathbf{Z}}_k^K{\mathbf{\Delta }}_k^K,{\mathbf{Z}}_k^K)$ using Algorithm 2

12:

compute [d_β]=LowRankDiag(XQ_β,X),

${\mathit{d}}_{\mathbf{Q}}=[\mathrm{diag}{\mathbf{Q}}_{\mathrm{s}}+{\mathit{\alpha }}^{-\mathrm{2}}{\mathit{d}}_{\mathit{\beta }};{\mathit{\alpha }}^{-\mathrm{2}}\mathrm{diag}{\mathbf{Q}}_{\mathit{\beta }}]$

13:

estimate posterior variance ${\mathit{d}}_{k,\mathit{\lambda }}={\mathit{\lambda }}^{-\mathrm{2}}{\mathit{d}}_{\mathbf{Q}}-{\mathit{d}}_{\mathrm{LR}}$

14:

return  MAP estimate γ_k,λ and variance estimate d_k,λ