Large-scale Inverse Problems Research Articles

Flexible Krylov Methods for Group Sparsity Regularization

Abstract This paper introduces new solvers for efficiently computing solutions to large-scale inverse problems with group sparsity regularization, including both non- overlapping and overlapping groups. Group sparsity regularization refers to a type of structured sparsity regularization, where the goal is to impose additional structure in the regularization process by assigning variables to predefined groups that may represent graph or network structures. Special cases of group sparsity regularization include l1 and isotropic total variation regularization. In this work, we develop hybrid projection methods based on flexible Krylov subspaces, where we first recast the group sparsity regularization term as a sequence of 2-norm penalization terms using adaptive regularization matrices in an iterative reweighted norm fashion. Then, we exploit flexible preconditioning techniques to efficiently incorporate the weight updates. The main advantages of these methods are that they are computationally efficient (leveraging the advantages of flexible methods), they are general (and therefore very easily adaptable to new regularization term choices), and they are able to select the regularization parameters automatically and adaptively (exploiting the advantages of hybrid methods). Extensions to multiple regularization terms and solution decomposition frameworks (e.g., for anomaly detection) are described, and a variety of numerical examples demonstrate both the efficiency and accuracy of the proposed approaches compared to existing solvers.

Efficient, multimodal, and derivative-free bayesian inference with Fisher–Rao gradient flows

Abstract In this paper, we study efficient approximate sampling for probability distributions known up to normalization constants. We specifically focus on a problem class arising in Bayesian inference for large-scale inverse problems in science and engineering applications. The computational challenges we address with the proposed methodology are: (i) the need for repeated evaluations of expensive forward models; (ii) the potential existence of multiple modes; and (iii) the fact that gradient of, or adjoint solver for, the forward model might not be feasible.

While existing Bayesian inference methods meet some of these challenges individually, we propose a framework that tackles all three systematically. Our approach builds upon the Fisher-Rao gradient flow in probability space, yielding a dynamical system for probability densities that converges towards the target distribution at a uniform exponential rate. This rapid convergence is advantageous for the computational burden outlined in (i). We apply Gaussian mixture approximations with operator splitting techniques to simulate the flow numerically; the resulting approximation can capture multiple modes thus addressing (ii). Furthermore, we employ the Kalman methodology to facilitate a derivative-free update of these Gaussian components and their respective weights, addressing the issue in (iii).

The proposed methodology results in an efficient derivative-free posterior approximation method, flexible enough to handle multi-modal distributions: Gaussian Mixture Kalman Inversion (GMKI). The effectiveness of GMKI is demonstrated both theoretically and numerically in several experiments with multimodal target distributions, including proof-of-concept and two-dimensional examples, as well as a large-scale application: recovering the Navier-Stokes initial condition from solution data at positive times.



A fast impact force identification method via constructing a dynamic reduced dictionary

Impact force identification is an important aspect of structural health monitoring for online assessment of the integrity of engineering structures. Increasing the number of monitoring positions is undoubtedly beneficial to accurate force localization, however, it may increase sharply the dimension of the transfer matrix and thus lead to a large-scale inverse problem. In this case, traditional methods, such as the Iterative Shrinkage Threshold algorithm (IST) and the Two-Step Iterative Shrinkage Thresholding algorithm (TwIST), suffer from low solving efficiency. In this paper, the dynamic reduced dictionary is firstly introduced according to the sparsity of impact force in both time and spatial domains so as to lower significantly the dimension of the transfer matrix. Then, a new reduced unconstrained optimization problem is established via the dynamic reduced dictionary for impact force identification. By integrating the reduced unconstrained optimization problem into each iteration step of IST/TwIST, two novel algorithms, which are the reduced Newton iteration shrinkage threshold algorithm (rNIST) and the reduced Newton two-step iteration shrinkage threshold algorithm (rNTwIST), are proposed to solve quickly the large-scale impact force identification problem using ℓ1-norm sparse regularization. Meanwhile, the adaptive setting of regularization parameters strategy is used to construct force more fast and accurately. The proposed algorithms are numerically and experimentally demonstrated through identifying impact forces at 16-point and 81-point distributions of a plate using one sensor response. The results show that rNIST/rNTwIST can more quickly and accurately identify impact forces under numerous monitoring points compared with other ℓ1-norm regularization methods and nonconvex sparse regularization methods. Moreover, the more the number of monitoring points, the more obvious the efficiency improvement of the impact force identification by rNIST/rNTwIST method.

TRIPs-Py: Techniques for regularization of inverse problems in python

AbstractIn this paper we describe TRIPs-Py, a new Python package of linear discrete inverse problems solvers and test problems. The goal of the package is two-fold: 1) to provide tools for solving small and large-scale inverse problems, and 2) to introduce test problems arising from a wide range of applications. The solvers available in TRIPs-Py include direct regularization methods (such as truncated singular value decomposition and Tikhonov) and iterative regularization techniques (such as Krylov subspace methods and recent solvers for $$\ell _p$$ ℓ p -$$\ell _q$$ ℓ q formulations, which enforce sparse or edge-preserving solutions and handle different noise types). All our solvers have built-in strategies to define the regularization parameter(s). Some of the test problems in TRIPs-Py arise from simulated image deblurring and computerized tomography, while other test problems model real problems in dynamic computerized tomography. Numerical examples are included to illustrate the usage as well as the performance of the described methods on the provided test problems. To the best of our knowledge, TRIPs-Py is the first Python software package of this kind, which may serve both research and didactical purposes.

Optimal design of large-scale nonlinear Bayesian inverse problems under model uncertainty

We consider optimal experimental design (OED) for Bayesian nonlinear inverse problems governed by partial differential equations (PDEs) under model uncertainty. Specifically, we consider inverse problems in which, in addition to the inversion parameters, the governing PDEs include secondary uncertain parameters. We focus on problems with infinite-dimensional inversion and secondary parameters and present a scalable computational framework for optimal design of such problems. The proposed approach enables Bayesian inversion and OED under uncertainty within a unified framework. We build on the Bayesian approximation error (BAE) approach, to incorporate modeling uncertainties in the Bayesian inverse problem, and methods for A-optimal design of infinite-dimensional Bayesian nonlinear inverse problems. Specifically, a Gaussian approximation to the posterior at the maximum a posteriori probability point is used to define an uncertainty aware OED objective that is tractable to evaluate and optimize. In particular, the OED objective can be computed at a cost, in the number of PDE solves, that does not grow with the dimension of the discretized inversion and secondary parameters. The OED problem is formulated as a binary bilevel PDE constrained optimization problem and a greedy algorithm, which provides a pragmatic approach, is used to find optimal designs. We demonstrate the effectiveness of the proposed approach for a model inverse problem governed by an elliptic PDE on a three-dimensional domain. Our computational results also highlight the pitfalls of ignoring modeling uncertainties in the OED and/or inference stages.

Uncertainty quantification for goal-oriented inverse problems via variational encoder-decoder networks

In this work, we describe a new approach that uses variational encoder-decoder (VED) networks for efficient uncertainty quantification for goal-oriented inverse problems. Contrary to standard inverse problems, these approaches are goal-oriented in that the goal is to estimate some quantities of interest (QoI) that are functions of the solution of an inverse problem, rather than the solution itself. Moreover, we are interested in computing uncertainty metrics associated with the QoI, thus utilizing a Bayesian approach for inverse problems that incorporates the prediction operator and techniques for exploring the posterior. This may be particularly challenging, especially for nonlinear, possibly unknown, operators and nonstandard prior assumptions. We harness recent advances in machine learning, i.e. VED networks, to describe a data-driven approach to large-scale inverse problems. This enables a real-time uncertainty quantification for the QoI. One of the advantages of our approach is that we avoid the need to solve challenging inversion problems by training a network to approximate the mapping from observations to QoI. Another main benefit is that we enable uncertainty quantification for the QoI by leveraging probability distributions in the latent and target spaces. This allows us to efficiently generate QoI samples and circumvent complicated or even unknown forward models and prediction operators. Numerical results from medical tomography reconstruction and nonlinear hydraulic tomography demonstrate the potential and broad applicability of the approach.

Parallel ensemble Kalman method with total variation regularization for large-scale field inversion

Field inversion is often encountered in data-driven computational modeling to infer latent spatial–varying parameters from available observations. The ensemble Kalman method is emerging as a useful tool for solving field inversion problems due to its derivative-free merits. However, the method is computationally prohibitive for large-scale field inversion with high-dimensional observation data, which necessitates developing a practical efficient implementation strategy. In this work, we propose a parallel implementation of the ensemble Kalman method with total variation regularization for large-scale field inversion problems. It is achieved by partitioning the computational domain into non-overlapping subdomains and performing local ensemble Kalman updates at each subdomain parallelly. In doing so, the computational complexity of the ensemble-based inversion method is significantly reduced to the level of local subdomains. Further, the total variation regularization is employed to smoothen the physical field over the entire domain, which can reduce the inference discrepancy caused by missing covariances near subdomain interfaces. The capability of the proposed method is demonstrated in three field inversion problems with increasing complexity, i.e., the diffusion problem, the scalar transport problem and the Reynolds averaged Navier-Stokes closure problem. The numerical results show that the proposed method can significantly improve computational efficiency with satisfactory inference accuracy.

Multilevel dimension-independent likelihood-informed MCMC for large-scale inverse problems

We present a non-trivial integration of dimension-independent likelihood-informed (DILI) MCMC (Cui et al 2016) and the multilevel MCMC (Dodwell et al 2015) to explore the hierarchy of posterior distributions. This integration offers several advantages: First, DILI-MCMC employs an intrinsic likelihood-informed subspace (LIS) (Cui et al 2014)—which involves a number of forward and adjoint model simulations—to design accelerated operator-weighted proposals. By exploiting the multilevel structure of the discretised parameters and discretised forward models, we design a Rayleigh–Ritz procedure to significantly reduce the computational effort in building the LIS and operating with DILI proposals. Second, the resulting DILI-MCMC can drastically improve the sampling efficiency of MCMC at each level, and hence reduce the integration error of the multilevel algorithm for fixed CPU time. Numerical results confirm the improved computational efficiency of the multilevel DILI approach.

A dimension-reduced variational approach for solving physics-based inverse problems using generative adversarial network priors and normalizing flows

We propose a novel modular inference approach combining two different generative models — generative adversarial networks (GAN) and normalizing flows — to approximate the posterior distribution of physics-based Bayesian inverse problems framed in high-dimensional ambient spaces. We dub the proposed framework GAN-Flow. The proposed method leverages the intrinsic dimension reduction and superior sample generation capabilities of GANs to define a low-dimensional data-driven prior distribution. Once a trained GAN-prior is available, the inverse problem is solved entirely in the latent space of the GAN using variational Bayesian inference with normalizing flow-based variational distribution, which approximates low-dimensional posterior distribution by transforming realizations from the low-dimensional latent prior (Gaussian) to corresponding realizations of a low-dimensional variational posterior distribution. The trained GAN generator then maps realizations from this approximate posterior distribution in the latent space back to the high-dimensional ambient space. We also propose a two-stage training strategy for GAN-Flow wherein we train the two generative models sequentially. Thereafter, GAN-Flow can estimate the statistics of posterior-predictive quantities of interest at virtually no additional computational cost. The synergy between the two types of generative models allows us to overcome many challenges associated with the application of Bayesian inference to large-scale inverse problems, chief among which are describing an informative prior and sampling from the high-dimensional posterior. GAN-Flow does not involve Markov chain Monte Carlo simulation, making it particularly suitable for solving large-scale inverse problems. We demonstrate the efficacy and flexibility of GAN-Flow on various physics-based inverse problems of varying ambient dimensionality and prior knowledge using different types of GANs and normalizing flows. Notably, one of the applications we consider involves a 65,536-dimensional inverse problem of phase retrieval wherein an object is reconstructed from sparse noisy measurements of the magnitude of its Fourier transform.

An optimized, easy-to-use, open-source GPU solver for large-scale inverse homogenization problems

An optimized, easy-to-use, open-source GPU solver for large-scale inverse homogenization problems

Hybrid Projection Methods for Solution Decomposition in Large-Scale Bayesian Inverse Problems

Hybrid Projection Methods for Solution Decomposition in Large-Scale Bayesian Inverse Problems

A differentiable, physics-informed ecosystem modeling and learning framework for large-scale inverse problems: demonstration with photosynthesis simulations

Abstract. Photosynthesis plays an important role in carbon, nitrogen, and water cycles. Ecosystem models for photosynthesis are characterized by many parameters that are obtained from limited in situ measurements and applied to the same plant types. Previous site-by-site calibration approaches could not leverage big data and faced issues like overfitting or parameter non-uniqueness. Here we developed an end-to-end programmatically differentiable (meaning gradients of outputs to variables used in the model can be obtained efficiently and accurately) version of the photosynthesis process representation within the Functionally Assembled Terrestrial Ecosystem Simulator (FATES) model. As a genre of physics-informed machine learning (ML), differentiable models couple physics-based formulations to neural networks (NNs) that learn parameterizations (and potentially processes) from observations, here photosynthesis rates. We first demonstrated that the framework was able to correctly recover multiple assumed parameter values concurrently using synthetic training data. Then, using a real-world dataset consisting of many different plant functional types (PFTs), we learned parameters that performed substantially better and greatly reduced biases compared to literature values. Further, the framework allowed us to gain insights at a large scale. Our results showed that the carboxylation rate at 25 ∘C (Vc,max25) was more impactful than a factor representing water limitation, although tuning both was helpful in addressing biases with the default values. This framework could potentially enable substantial improvement in our capability to learn parameters and reduce biases for ecosystem modeling at large scales.

A variable projection method for large-scale inverse problems with ℓ1 regularization

Inference by means of mathematical modeling from a collection of observations remains a crucial tool for scientific discovery and is ubiquitous in application areas such as signal compression, imaging restoration, and supervised machine learning. With ever-increasing model complexities and larger data sets, new specially designed methods are urgently needed to recover meaningful quantities of interest. We consider the broad spectrum of linear inverse problems where the aim is to reconstruct quantities with a sparse representation on some vector space. We provide a new variable projection augmented Lagrangian algorithm to solve the underlying ℓ1 regularized inverse problem that is both efficient and effective. We present the proof of convergence for an algorithm using an inexact step for the projected problem at each iteration. The performance and convergence properties for various imaging problems are investigated. The efficiency of the algorithm makes it feasible to automatically find the regularization parameter, here illustrated using an argument based on the degrees of freedom of the objective function equipped with a bisection algorithm for root-finding.

Hierarchical off-diagonal low-rank approximation of Hessians in inverse problems, with application to ice sheet model initialization

Obtaining lightweight and accurate approximations of discretized objective functional Hessians in inverse problems governed by partial differential equations (PDEs) is essential to make both deterministic and Bayesian statistical large-scale inverse problems computationally tractable. The cubic computational complexity of dense linear algebraic tasks, such as Cholesky factorization, that provide a means to sample Gaussian distributions and determine solutions of Newton linear systems is a computational bottleneck at large-scale. These tasks can be reduced to log-linear complexity by utilizing hierarchical off-diagonal low-rank (HODLR) matrix approximations. In this work, we show that a class of Hessians that arise from inverse problems governed by PDEs are well approximated by the HODLR matrix format. In particular, we study inverse problems governed by PDEs that model the instantaneous viscous flow of ice sheets. In these problems, we seek a spatially distributed basal sliding parameter field such that the flow predicted by the ice sheet model is consistent with ice sheet surface velocity observations. We demonstrate the use of HODLR Hessian approximation to efficiently sample the Laplace approximation of the posterior distribution with covariance further approximated by HODLR matrix compression. Computational studies are performed which illustrate ice sheet problem regimes for which the Gauss–Newton data-misfit Hessian is more efficiently approximated by the HODLR matrix format than the low-rank (LR) format. We then demonstrate that HODLR approximations can be favorable, when compared to global LR approximations, for large-scale problems by studying the data-misfit Hessian associated with inverse problems governed by the first-order Stokes flow model on the Humboldt glacier and Greenland ice sheet.

A differentiable, physics-informed ecosystem modeling and learning framework for large-scale inverse problems

A differentiable, physics-informed ecosystem modeling and learning framework for large-scale inverse problems

HIPPYlib-MUQ: A Bayesian Inference Software Framework for Integration of Data with Complex Predictive Models under Uncertainty

Bayesian inference provides a systematic framework for integration of data with mathematical models to quantify the uncertainty in the solution of the inverse problem. However, the solution of Bayesian inverse problems governed by complex forward models described bypartial differential equations (PDEs)remains prohibitive with black-boxMarkov chain Monte Carlo (MCMC)methods. We present hIPPYlib-MUQ, an extensible and scalable software framework that contains implementations of state-of-the art algorithms aimed to overcome the challenges of high-dimensional, PDE-constrained Bayesian inverse problems. These algorithms accelerate MCMC sampling by exploiting the geometry and intrinsic low-dimensionality of parameter space via derivative information and low rank approximation. The software integrates two complementary open-source software packages, hIPPYlib and MUQ. hIPPYlib solves PDE-constrained inverse problems using automatically-generated adjoint-based derivatives, but it lacks full Bayesian capabilities. MUQ provides a spectrum of powerful Bayesian inversion models and algorithms, but expects forward models to come equipped with gradients and Hessians to permit large-scale solution. By combining these two complementary libraries, we created a robust, scalable, and efficient software framework that realizes the benefits of each and allows us to tackle complex large-scale Bayesian inverse problems across a broad spectrum of scientific and engineering disciplines. To illustrate the capabilities of hIPPYlib-MUQ, we present a comparison of a number of MCMC methods available in the integrated software on several high-dimensional Bayesian inverse problems. These include problems characterized by both linear and nonlinear PDEs, various noise models, and different parameter dimensions. The results demonstrate that large (∼ 50×) speedups over conventional black box and gradient-based MCMC algorithms can be obtained by exploiting Hessian information (from the log-posterior), underscoring the power of the integrated hIPPYlib-MUQ framework.

A differentiable, physics-informed ecosystem modeling and learning framework for large-scale inverse problems

A differentiable, physics-informed ecosystem modeling and learning framework for large-scale inverse problems

Overlapping Domain Decomposition Method with Cascadic Multigrid for Image Restoration

In the process of image restoration, it is usually necessary to solve large-scale inverse problems, where the computational cost is very high for large or high-resolution images. The domain decomposition method is one of the most effective algorithms to solve large-scale problems, which can effectively decrease the computational cost. The cascadic multigrid method has a good effect on the linear model of image restoration and can obtain high quality restored images. In this paper, the overlapping domain decomposition method (DDM) with the cascadic multigrid method (CMG) and the DDM with new extrapolation cascadic multigrid method (NECMG) are presented to solve the image restoration problems of denoising and deblurring. We first divide the image problem into some overlapping and independent subproblems. Then, each subproblem is solved independently by CMG or NECMG with the edge-preserving operator. Numerical experiments show that the new method is effective.

On global normal linear approximations for nonlinear Bayesian inverse problems

The replacement of a nonlinear parameter-to-observable mapping with a linear (affine) approximation is often carried out to reduce the computational costs associated with solving large-scale inverse problems governed by partial differential equations (PDEs). In the case of a linear parameter-to-observable mapping with normally distributed additive noise and a Gaussian prior measure on the parameters, the posterior is Gaussian. However, substituting an accurate model for a (possibly well justified) linear surrogate model can give misleading results if the induced model approximation error is not accounted for. To account for the errors, the Bayesian approximation error (BAE) approach can be utilised, in which the first and second order statistics of the errors are computed via sampling. The most common linear approximation is carried out via linear Taylor expansion, which requires the computation of (Fréchet) derivatives of the parameter-to-observable mapping with respect to the parameters of interest. In this paper, we prove that the (approximate) posterior measure obtained by replacing the nonlinear parameter-to-observable mapping with a linear approximation is in fact independent of the choice of the linear approximation when the BAE approach is employed. Thus, somewhat non-intuitively, employing the zero-model as the linear approximation gives the same approximate posterior as any other choice of linear approximations of the parameter-to-observable model. The independence of the linear approximation is demonstrated mathematically and illustrated with two numerical PDE-based problems: an inverse scattering type problem and an inverse conductivity type problem.

Uncertainty Quantification and Experimental Design for Large-Scale Linear Inverse Problems under Gaussian Process Priors

We consider the use of Gaussian process (GP) priors for solving inverse problems in a Bayesian framework. As is well known, the computational complexity of GPs scales cubically in the number of datapoints. Here we show that in the context of inverse problems involving integral operators, one faces additional difficulties that hinder inversion on large grids. Furthermore, in that context, covariance matrices can become too large to be stored. By leveraging recent results about sequential disintegrations of Gaussian measures, we are able to introduce an implicit representation of posterior covariance matrices that reduces the memory footprint by only storing low rank intermediate matrices, while allowing individual elements to be accessed on-the-fly without needing to build full posterior covariance matrices. Moreover, it allows for fast sequential inclusion of new observations. These features are crucial when considering sequential experimental design tasks. We demonstrate our approach by computing sequential data collection plans for excursion set recovery for a gravimetric inverse problem, where the goal is to provide fine resolution estimates of high density regions inside the Stromboli volcano, Italy. Sequential data collection plans are computed by extending the weighted integrated variance reduction (wIVR) criterion to inverse problems. Our results show that this criterion is able to significantly reduce the uncertainty on the excursion volume, reaching close to minimal levels of residual uncertainty. Overall, our techniques allow the advantages of probabilistic models to be brought to bear on large-scale inverse problems arising in the natural sciences. Particularly, applying the latest developments in Bayesian sequential experimental design on realistic large-scale problems opens new venues of research at a crossroads between mathematical modelling of natural phenomena, statistical data science, and active learning.