4.1 Quasi-Newton approach

Many existing variational or adjoint-based inverse modeling studies use an iterative, quasi-Newton approach to estimate the fluxes (e.g., Baker et al., 2006; Henze et al., 2007). This strategy requires creating scripts that calculate the cost function (Eq. 6) and gradient, and the quasi-Newton solver will iterate toward the solution using these two inputs. The most common quasi-Newton solver in existing studies is the limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS) (Nocedal, 1980; Liu and Nocedal, 1989). Furthermore, one variant of this algorithm, L-BFGS-B, will estimate the fluxes (s) subject to a bound. For example, some trace gases do not have large surface sinks, and L-BFGS-B can ensure that the estimated fluxes are non-negative. Refer to Miller et al. (2014) for a full discussion of strategies to enforce bounds on atmospheric inverse problems.

The implementation of L-BFGS or L-BFGS-B for a GIM is more complicated than in a classical Bayesian inverse model. Specifically, the goal is to use one of these algorithms to estimate the fluxes (s), but the cost function in the GIM has an additional unknown variable (β, Eq. 6). We substitute the equation for the unknown coefficients (β) into the GIM cost function (Eq. 6), thereby removing β (e.g., Kitanidis, 1995):

where G has dimensions m×m. Note that we never formulate G explicitly when calculating Eq. (11) because it is often a large, non-sparse m×m matrix. Rather, we successively multiply s by the individual components of Eq. (12) to avoid formulating a full m×m matrix. In addition, we use the Kronecker product identity in Eq. (10) to avoid explicitly formulating Q⁻¹. The SI describes this approach in greater detail.

We can further speed up the convergence of L-BFGS using a variable transformation. We transform the fluxes from s to s^*, similar to the approach used in a handful of atmospheric inverse modeling studies (e.g., Baker et al., 2006; Meirink et al., 2008). The strategy is to first solve a transformed optimization problem for s^* and then to obtain s using the following relations:

where Q^½ is the symmetric square root of Q. Note that we never explicitly formulate Q^½ but instead do all matrix operations on the individual components of Q^½. We then substitute the equation for s^* into the cost function (L_s):

The functions $L_{{\mathit{s}}^{*}}$ and ∇L(s^*) are then used as inputs to the L-BFGS algorithm, and the resulting optimization can converge in fewer iterations than using the cost function without the transformation (Eq. 11). Section 7 includes further discussion of why L-BFGS without a transformation may converge slowly and how a variable transformation can remedy the problem.

4.2 Minimum-residual approach

The minimum-residual approach described here uses a very different strategy for the iterative optimization (Paige and Saunders, 1975). The quasi-Newton approach described above will search for the minimum of the cost function with the help of the gradient (Eqs. 16–17). By contrast, there is a class of solvers that will estimate the solution to a linear system of equations where one side is too large to invert, and this class of solvers offers an alternative strategy. Specifically, these methods will solve a system of equations of the form Ax=b, where x is a vector of unknown values, b is a known vector, and A is a matrix that is too large to store in memory or too expensive to form explicitly. This strategy has been employed by inverse modeling studies in hydrology and seismology (Saibaba and Kitanidis, 2012; Saibaba et al., 2012; Liu et al., 2014b; Lee et al., 2016).

This strategy can also be employed to estimate the fluxes using Eq. (8). We cannot solve these equations directly for inverse problems with large datasets (i.e., large n); the left-hand side of Eq. (8) becomes too large to invert and too expensive to explicitly form. Instead, we use the minimum-residual method to estimate ξ and β (e.g., Barrett et al., 1994). This class of algorithms only require a function that will calculate the left-hand side of Eq. (8) given some guess for ξ and β. Then, one can iteratively compute an approximate solution by forming a series of matrix–vector products – the product of ξ and β with the matrices on the left-hand side of Eq. (8). This approach also makes it feasible to estimate the fluxes (s) even if the atmospheric model does not explicitly calculate H or H^T. One can pass the vector ξ through the model adjoint to calculate H^Tξ and then pass the vector QH^Tξ through the forward model to calculate H(QH^Tξ). As a result, one never needs to calculate or store the $(n+p)\times (n+p)$ matrix on the left-hand side of Eq. (8) or invert that matrix.

Note that some studies use a preconditioner to help the minimum-residual algorithm converge more quickly to a solution (e.g., Saibaba and Kitanidis, 2012; Liu et al., 2014b; Lee et al., 2016). The preconditioner can speed up convergence by reducing the condition number and/or clustering the eigenvalues near 1. Saibaba and Kitanidis (2012) detail one possible strategy for preconditioning the GIM, and this approach has subsequently been employed in several studies (e.g., Liu et al., 2014b; Lee et al., 2016). We do not detail the implementation here, but Saibaba and Kitanidis (2012) describe the step-by-step procedure. The preconditioner detailed in that study can dramatically speed up convergence (Sect. 7) but may require tens to thousands of forward model runs as a precomputational cost to realize these improvements in convergence, depending upon the complexity of the covariance matrix Q. Note that these forward model runs can be done simultaneously in parallel. A discussion of the computational cost is provided in Sect. 7.

5.1 Conditional realizations or simulations

One approach to estimate the posterior uncertainties is to generate conditional realizations or Monte Carlo simulations, and several variational or adjoint-based studies of greenhouse gases have employed this strategy (e.g., Chevallier et al., 2007; Liu et al., 2014a; Bousserez et al., 2015). A conditional realization is a random sample from the posterior distribution (e.g., Kitanidis, 1995; Michalak et al., 2004; Chevallier et al., 2007), and the statistics of many conditional realizations will approximate the posterior variances and covariances. These variances and covariances can have complex spatial and temporal patterns, and it can be challenging to adequately sample the tails of the posterior distribution. As a result, several studies of atmospheric trace gases use thousands of realizations to sample the posterior distribution and approximate the uncertainties (e.g., Rigby et al., 2011; Ganesan et al., 2014). In other cases, computational constraints make it impractical to generate more than tens or hundreds of realizations, especially for very large variational or adjoint-based inverse problems (e.g., Chevallier et al., 2007; Liu et al., 2014a).

One approach for generating a single conditional realization is similar to the procedure for calculating the best estimate of the fluxes. To generate a conditional realization, we first generate random samples from 𝒩(0,R) and 𝒩(0,Q) and use these samples to perturb the data and the parameters respectively; we then solve the GIM using any of the three approaches described previously. In each case, the equations are slightly different, and we therefore do not list all of the equations here. Michalak et al. (2004) and Saibaba and Kitanidis (2012) describe how to generate conditional realizations using the direct approach in Sect. 3, Kitanidis (1995) and Snodgrass and Kitanidis (1997) and the SI describe how to generate realizations using quasi-Newton methods (Sect. 4.1), and Kitanidis (1996) and Saibaba and Kitanidis (2012) describe how to generate realizations using the minimum-residual approach (Sect. 4.2).

5.2 Reduced-rank approach

A number of studies reduce the computational burden of uncertainty estimation by replacing one matrix in the calculations with a low-rank approximation. This strategy is based upon the direct calculation of the posterior covariance matrix (${\mathbf{V}}_{\widehat{\mathit{s}}}$) (e.g., Saibaba and Kitanidis, 2015):

where ${\mathbf{V}}_{\widehat{\mathit{s}}}$ and V₁ have dimensions m×m, V₂ dimensions m×p, and V₃ dimensions p×p. Note that V₁ is the posterior covariance matrix in a classical Bayesian inverse model. The uncertainty calculations in a GIM include an additional term, notably ${\mathbf{V}}_{\mathrm{2}}{\mathbf{V}}_{\mathrm{3}}{\mathbf{V}}_{\mathrm{2}}^T$. This term accounts for the effect of uncertain coefficients (β) on the estimated fluxes. These coefficients are estimated as part of the GIM, and uncertainty in these coefficients contributes additional uncertainty to the posterior flux estimate ($\widehat{\mathit{s}}$). Also note that there are many equivalent equations for calculating the posterior covariance matrix, and the form in Eq. (19) is particularly conducive to a low-rank approximation strategy.

This direct calculation of the posterior covariance matrix is not tractable for large inverse problems, in part, because it requires inverting the matrix sum in V₁. Furthermore, it requires computing the matrix–matrix product H^TR⁻¹H, a step that can be impractical for models that do not explicitly formulate H.

Instead, we use a specific low-rank approximation in the calculations for V₁. In other words, we approximate a matrix using a limited number of eigenvectors and eigenvalues, thereby making the calculations for V₁ computationally feasible. A number of studies discuss this strategy in the context of classical Bayesian inverse modeling (e.g., Meirink et al., 2008; Flath et al., 2011; Spantini et al., 2015). Here, we review how this approach can be applied in the context of a GIM (e.g., Saibaba and Kitanidis, 2015). We summarize the procedure here, but refer to Flath et al. (2011), Spantini et al. (2015), and Saibaba and Kitanidis (2015) for additional detail. The main idea is to consider the matrix

where Q^½ can be the symmetric square root or Cholesky decomposition of Q. This matrix is sometimes referred to as the prior-preconditioned data-misfit part of the Hessian. Previous studies have leveraged the fact that, in many applications, this matrix has rapidly decaying eigenvalues and therefore, can be accurately approximated using a low-rank matrix (e.g., Meirink et al., 2008; Flath et al., 2011; Spantini et al., 2015). This low-rank representation has a dual purpose: it is more efficient to store this matrix in memory and improves the efficiency of computations. There are several methods to compute this low-rank approximation, and these methods can estimate the eigenvalues and vectors without storing the above matrix (Eq. 23) in memory. Instead, these algorithms (e.g., Krylov subspace methods) iteratively estimate the largest eigenvalues and eigenvectors of a matrix A given some function that calculates Ax, where x is a vector that is provided by the algorithm and does not need to be specified by the user. For example, the eigs function in MATLAB offers an interface to these algorithms.

In this study, we estimate the eigenvectors and eigenvalues using a randomized algorithm developed in Halko et al. (2011). This approach requires just over two forward model runs and two adjoint model runs per (approximate) eigenpair, and this algorithm is therefore less computationally intensive than many other available algorithms. Furthermore, these forward and adjoint model runs can be generated in parallel, reducing the required computing time. Using this approach, we obtain the low-rank approximation

where U (n×ℓ) contains the approximate eigenvectors, Λ (ℓ×ℓ) is a diagonal matrix whose diagonals contain the approximate eigenvalues, and ℓ is the number of approximate eigenpairs computed. Note that Saibaba and Kitanidis (2015) and Saibaba et al. (2016) provide an alternative approach that does not require taking the symmetric square root or Cholesky decomposition of Q. That approach is a good choice if the symmetric square root or Cholesky decomposition of D and/or E are difficult to compute.

We can use these eigenvalues or vectors to approximate V₁ using the Woodbury matrix identity (Flath et al., 2011; Spantini et al., 2015):

In this setup, V₁ is written as an update to the prior covariance matrix. An intuitive way of understanding is that by observing data, the variance (i.e., the uncertainty) is reduced since we know more about the fluxes or emissions of interest. The prior covariance matrix (Q) in Eq. (25) is taken to be a positive definite, full-rank matrix and is not affected at all by the approximation. Rather, the update is approximated using a limited number of eigenvectors and eigenvalues. We subtract this positive semidefinite update, ensuring that the variance in V₁ is reduced. The more eigenvectors and eigenvalues used in Eq. (25), the better the approximation of this update. The estimate for V₁ can subsequently be plugged into the expression for ${\mathbf{V}}_{\widehat{\mathit{s}}}$ to obtain an approximation for the posterior covariance matrix.

This reduced-rank approach greatly improves the tractability of posterior uncertainty calculations, but one computational road block remains; ${\mathbf{V}}_{\widehat{\mathit{s}}}$ is usually too large to store in memory. However, one can calculate uncertainties for individual grid boxes or for aggregate regions without storing ${\mathbf{V}}_{\widehat{\mathit{s}}}$. To do so, multiply the right-hand side of Eq. (19) by a vector, resulting in a series of matrix–vector calculations. In this case, the vector should have a 1 in the flux grid box(es) of interest and a 0 elsewhere. If that vector has a 1 for multiple elements, then the calculation will yield the uncertainty in the flux estimate summed across several model grid boxes. For example, the following equation will compute the uncertainty in the total flux summed over all locations and times:

where 1 (m×1) is a vector of ones. Note that it may be necessary to convert the units of the fluxes in the course of these calculations. For example, atmospheric models that use a latitude–longitude grid will have grid boxes with different areas, and it may be necessary to multiply the vector 1 by the area of each grid box to account for these differing grid box areas.

7.1 Best estimate of the fluxes

We test several algorithms for estimating the fluxes (s), and all but one of these algorithms converges quickly toward the solution. The minimum-residual approach (with and without a preconditioner) and the L-BFGS with a variable transformation converge quickly for both case studies (Figs. 1, 2, and 3). Note that we cannot compare the larger case study to a direct solution, but we do compare the results against CO₂ fluxes estimated using a very large number of iterations (250 in this case).

Figure 1Root mean squared error (RMSE) relative to the direct solution for the 6-week case study. The top panel (a) displays RMSE calculated using each of the 1.05×10⁶ flux grid boxes. The minimum-residual method with a preconditioner converges most quickly on the direct solution while L-BFGS shows poor convergence. We also average the 3-hourly posterior flux estimate to obtain a monthly averaged flux for each model grid box, and the error in this monthly average is displayed in panel (b). The flux estimate, when averaged up to this aggregate monthly scale, converges more quickly to the direct solution than the individual 3-hourly fluxes.

Download

Figure 2The “true” CO₂ fluxes from CT2017 used in the 6-week case study (a) and the direct estimate of the fluxes (b). We estimate 3-hourly fluxes but average them across the month of July in this figure. Panels (c–e) display the flux estimate from L-BFGS as a function of iteration number, panels (f–h) those from L-BFGS with a variable transformation, and panels (i–k) those from the minimum-residual approach. Both L-BFGS with the transformation and the minimum-residual approach resemble the direct solution after about 20 iterations (panels d, g, and j).

Figure 3RMSE relative to the direct solution for the full-year case study. The top panel (a) displays RMSE calculated using each of the 9.4×10⁶ fluxes. We also average the 3-hourly posterior flux estimate to obtain an annually averaged flux for each model grid box (b). The full-year case study is much larger than the 6-week case study (i.e., includes more observations and unknowns), yet the flux estimate converges to the solution after a similar number of iterations.

Download

By contrast, the L-BFGS algorithm without a variable transformation converges very slowly to the solution (Figs. 1 and 3). The L-BFGS algorithm makes a simple approximation for the posterior covariance matrix at each iteration (i.e., the inverse Hessian); it uses this approximation and the gradient of the cost function to determine the direction of steepest descent and iterate toward the solution (e.g., Nocedal, 1980; Liu and Nocedal, 1989). Most L-BFGS algorithms use a memory-saving sparse matrix like the identity matrix as an initial guess for the inverse Hessian. In the case studies here, the inverse Hessian without the variable transformation has complex off-diagonal elements, and the identity matrix is therefore a poor approximation. By contrast, the inverse Hessian with the variable transformation has relatively small off-diagonal elements, and the identity matrix is a much better approximation. This difference likely explains why the L-BFGS algorithm without the variable transformation converges far more slowly than the algorithm with the transformation.

Of all algorithms, the minimum-residual algorithm with preconditioning converges in the fewest number of iterations (Fig. 1). However, this faster convergence comes at a cost. We test a preconditioner with ℓ=400 and ℓ=2500 approximate eigenvectors. Only the preconditioner with 2500 eigenvectors yields faster convergence and only in the smaller of the two case studies. To construct this preconditioner, we had to apply each of the 2500 eigenvectors to the forward model (H). These forward model runs may be achievable if an explicit H matrix is available, if the forward model is not computationally expensive, and/or if a large computing cluster is available to distribute these runs across many cores. However, there are many instances when it may be impractical to generate such a large number of forward model runs. Since the preconditioner can be used for computing the best estimate as well as generating the conditional realizations, the upfront cost of constructing the preconditioner may be reutilized in the uncertainty quantification step in the following ways: (1) the cost can be amortized if hundreds or thousands of conditional realizations are generated, and (2) information collected during the construction can also be used in the reduced-rank computations. Moreover, the use of the randomized algorithm for constructing the preconditioners can be parallelized across the forward model runs. It may also be possible to develop a more computationally efficient preconditioner, but that objective is beyond the scope of this study.

It is also important to note that these algorithms converge more quickly on an accurate flux estimate for aggregate space and timescales. In this study, we always estimate the fluxes at a 3-hourly time resolution. After estimating the fluxes, we average the estimate within each grid box across the entire month or year and compare this time-averaged flux to the direct solution (Figs. 1 and 3). These monthly and annual averages are comparable to the direct solution after only a few iterations of the L-BFGS or minimum-residual algorithms. By contrast, it takes more iterations for each algorithm to converge on a flux estimate that is accurate relative to the direct solution for each 3 h time interval. Many GIM studies of CO₂ report monthly flux totals (e.g., Gourdji et al., 2012; Shiga et al., 2018), and these monthly totals may therefore be a more important quantity to robustly estimate than 3-hourly fluxes.

Furthermore, the number of iterations required to converge on the solution does not appear to change dramatically with the size of the problem (Figs. 1 and 3). Most algorithms converge for aggregate time periods (Fig. 1b) after 50 iterations in the 6-week case study and after about 75 iterations in the year-long case study (Fig. 3b).

Ultimately, the best or optimal algorithm will likely depend upon the specifics of the inverse modeling problem in question. For example, several existing variational or adjoint-based inverse models are already built to work with the L-BFGS algorithm, and it may be more convenient to adapt this existing infrastructure to the geostatistical approach (e.g., Baker et al., 2006; Henze et al., 2007). By contrast, the minimum-residual approach does not require inverting or decomposing the covariance matrix components (e.g., R, D, E), and this approach is advantageous if any of these components is difficult to invert. For example, existing inverse modeling studies using OCO-2 observations construct R as a diagonal matrix, comparable to the setup here (e.g., Crowell et al., 2019). In future applications, it may be important to include off-diagonal elements or spatial and temporal error covariances in R – to better capture the information content of the satellite observations and to estimate robust uncertainty bounds, among other reasons. The quasi-Newton approach using L-BFGS requires inverting R, but this calculation may not be computationally feasible if R has off-diagonal elements. By contrast, the minimum-residual approach does not require the inverse of R. As a result, one could compute elements of R on the fly, without ever storing R in its entirety. In addition to these considerations, the choice of whether or not to use a preconditioner with minimum-residual approach depends upon the computational burden of passing many vectors through the forward model. However, the development of new preconditioning strategies could alter this trade-off.

7.2 Uncertainty quantification

The reduced-rank approach produces the most conservative uncertainty estimates; it will typically approximate uncertainties that are equal to or larger than the true uncertainties. By contrast, the conditional realizations will typically underestimate the posterior uncertainties. Figures 4 and 5 show the uncertainties estimated for both case studies as a function of the number of forward or adjoint model runs. Both approaches converge toward the true posterior uncertainties as the number of eigenvectors or conditional realizations increases.

Figure 4Estimated uncertainties using conditional realizations and the reduced-rank approach for the 6-week case study. Panel (a) shows the estimated uncertainty in the total estimated CO₂ flux for the US and panel (b) for New England. Both approaches converge on the direct estimate as the number of realizations and eigenpairs increase. Note that we estimate the uncertainties for each model grid box and each 3 h time period, but we report the uncertainties for broader aggregate regions in this figure. Most existing GIM studies report the uncertainties for large regions because these regions are either more ecologically or policy relevant than the individual model grid boxes. Additionally, the uncertainties are almost always smaller proportional to the total flux for large regions than for individual grid boxes; atmospheric observations usually provide a much more robust constraint on regional flux totals than for individual model grid boxes.

Download

Figure 5Estimated uncertainties for the full-year case study, analogous to Fig. 4. The uncertainties estimated by the reduced-rank approach are consistently higher than the conditional realizations, though they begin to converge as the number of realizations and eigenpairs increases. Note that it is not possible to obtain a direct estimate for the uncertainties given the size of this case study, so we cannot evaluate the accuracy of the uncertainty estimates shown here.

Download

These approaches will tend to under- and overestimate the uncertainties due to the approximations involved in each approach. Conditional realizations randomly sample from the posterior uncertainties, and one may need to generate hundreds or thousands of realizations to effectively sample the entire uncertainty space, particularly for large, complex inverse problems. By contrast, the reduced-rank approach approximates the posterior uncertainties by subtracting an update from the prior covariance matrix (Q) (Eq. 25). This update will always be too small if it is approximated using a limited number of eigenvectors and eigenvalues, yielding posterior uncertainties that are too large.

Neither of these approaches provides a silver bullet, so to speak, for estimating the uncertainties. Both require a large number of forward and adjoint model runs, a requirement that can be computationally intensive, depending on the requirements for the forward and adjoint models. These model runs can be generated in parallel for either approach, ameliorating this computational burden. By contrast, one can produce an uncertainty estimate using a smaller number of model runs, but the resulting uncertainty estimates will either be too small (using conditional realizations) or too large (using the reduced-rank approach). It is arguably more conservative to generate posterior uncertainties that are too large than too small; the latter may cause a modeler to overstate the results or incorrectly conclude that a result is statistically significant when it is not. However, neither of these outcomes is ideal.