High-dimensional Inverse Problems Research Articles

VIP - Variational Inversion Package with example implementations of Bayesian tomographic imaging

Bayesian inference has become an important methodology to solve inverse problems and to quantify uncertainties in their solutions. Variational inference is a method that provides probabilistic, Bayesian solutions efficiently by using optimisation. In this study we present a Python Variational Inversion Package (VIP), to solve inverse problems using variational inference methods. The package includes automatic differential variational inference (ADVI), Stein variational gradient descent (SVGD) and stochastic SVGD (sSVGD), and provides implementations of 2D travel time tomography and 2D full waveform inversion including test examples and solutions. Users can solve their own problems by supplying an appropriate forward function and a gradient calculation code. In addition, the package provides a scalable implementation which can be deployed easily on a desktop machine or using modern high performance computational facilities. The examples demonstrate that VIP is an efficient, scalable, extensible and user-friendly package, and can be used to solve a wide range of low or high dimensional inverse problems in practice.

Maximum-likelihood estimators in physics-informed neural networks for high-dimensional inverse problems

Physics-informed neural networks (PINNs) have proven a suitable mathematical scaffold for solving inverse ordinary (ODE) and partial differential equations (PDE). Typical inverse PINNs are formulated as soft-constrained multi-objective optimization problems with several hyperparameters. In this work, we demonstrate that inverse PINNs can be framed in terms of maximum-likelihood estimators (MLE) to allow explicit error propagation from interpolation to the physical model space through a Taylor expansion, removing a hyperparameter from the formulation. We demonstrate the technique through its application to high-dimensional coupled ODEs constrained by differential algebraic equations that are common in transient chemical and biological kinetics. Furthermore, we show that singular-value decomposition (SVD) of the ODE coupling matrix (reaction stoichiometry matrix) provides reduced uncorrelated subspaces in which PINNs solutions can be represented and over which residuals can be projected. Finally, SVD bases serve as preconditioners for the inversion of covariance matrices in this hyperparameter-free robust application of MLE to “kinetics-informed neural networks”.

A Distributed Block-Split Gibbs Sampler with Hypergraph Structure for High-Dimensional Inverse Problems

Sampling-based algorithms are classical approaches to perform Bayesian inference in inverse problems. They provide estimators with the associated credibility intervals to quantify the uncertainty on the estimators. Although these methods hardly scale to high dimensional problems, they have recently been paired with optimization techniques, such as proximal and splitting approaches, to address this issue. Such approaches pave the way to distributed samplers, splitting computations to make inference more scalable and faster. We introduce a distributed Split Gibbs sampler (SGS) to efficiently solve such problems involving distributions with multiple smooth and non-smooth functions composed with linear operators. The proposed approach leverages a recent approximate augmentation technique reminiscent of primal-dual optimization methods. It is further combined with a block-coordinate approach to split the primal and dual variables into blocks, leading to a distributed block-coordinate SGS. The resulting algorithm exploits the hypergraph structure of the involved linear operators to efficiently distribute the variables over multiple workers under controlled communication costs. It accommodates several distributed architectures, such as the Single Program Multiple Data and client-server architectures. Experiments on a large image deblurring problem show the performance of the proposed approach to produce high quality estimates with credibility intervals in a small amount of time. Supplementary material to reproduce the experiments is available online.

VI-DGP: A Variational Inference Method with Deep Generative Prior for Solving High-Dimensional Inverse Problems

VI-DGP: A Variational Inference Method with Deep Generative Prior for Solving High-Dimensional Inverse Problems

Comparative analyses of covariance matrix adaptation and iterative ensemble smoother on high-dimensional inverse problems in high-resolution groundwater modeling

Parameter estimation is an inverse problem which is crucial to reliable groundwater model predictions and management. Numerous techniques have been developed to address this challenging problem. This study aims to compare the performance of a stochastic optimization method, covariance matrix adaptation-evolution strategy (CMA-ES), and a data assimilation method, Levenberg-Marquardt based iterative ensemble smoother (ES-LM), on solving inverse problems in groundwater modeling. This study also presents a parallelization strategy to accelerate their implementation. The comparative analysis involves synthetic inverse problems and a real-world case study of calibrating a high-resolution groundwater model. Results from the synthetic problems suggest that both the CMA-ES and the ES-LM are able to achieve desirable data fitting, but the ES-LM generally exhibits greater efficiency than the CMA-ES. Also, with appropriate stopping criteria, both methods can perform well on propagating uncertainty of model predictions. However, both methods can underestimate the uncertainty of model parameters derived from Bayesian posterior. The underestimation is particularly pronounced by the ES-LM when using a small ensemble size, whereas uncertainty quantification by the CMA-ES is less affected by population size. Moreover, calibration of the highly parameterized groundwater model in Louisiana and Southwest Mississippi indicates that the ES-LM can outperform the CMA-ES for high-dimensional inverse problems regarding both data match and computational costs. With proper ensemble size, the two methods can produce comparable uncertainty of model parameters and predictions. This study contributes to the scientific understanding and efficient application of the CMA-ES and the ES-LM.

Machine learning based priors for Bayesian inversion in MR imaging

The Bayesian approach allows the incorporation of informative prior knowledge to effectively enable and improve the solution of inverse problems. Obtaining prior information in probabilistic terms is, however, a challenging task. Recently, machine learning has been applied for the training of generative models to facilitate the translation of historically or otherwise available data to a prior distribution. In this work, we apply this methodology to undersampled magnetic resonance imaging. In particular, we employ an autoencoder as part of a generative model to statistically regularise and solve the high-dimensional inverse problem using Bayesian inversion. Comparison with a classical Gaussian Markov random field prior is performed and numerical examples highlight the possible advantages of data-driven priors.

Learning to learn by using nonequilibrium training protocols for adaptable materials

Evolution in time-varying environments naturally leads to adaptable biological systems that can easily switch functionalities. Advances in the synthesis of environmentally responsive materials therefore open up the possibility of creating a wide range of synthetic materials which can also be trained for adaptability. We consider high-dimensional inverse problems for materials where any particular functionality can be realized by numerous equivalent choices of design parameters. By periodically switching targets in a given design algorithm, we can teach a material to perform incompatible functionalities with minimal changes in design parameters. We exhibit this learning strategy for adaptability in two simulated settings: elastic networks that are designed to switch deformation modes with minimal bond changes and heteropolymers whose folding pathway selections are controlled by a minimal set of monomer affinities. The resulting designs can reveal physical principles, such as nucleation-controlled folding, that enable such adaptability.

Localized ensemble Kalman inversion

Ensemble Kalman inversion (EKI) is an adaption of the ensemble Kalman filter (EnKF) for the numerical solution of inverse problems. Both EKI and EnKF suffer from the ‘subspace property’, i.e. the EKI and EnKF solutions are linear combinations of the initial ensembles. The subspace property implies that the ensemble size should be larger than the problem dimension to ensure EKI’s convergence to the correct solution. This scaling of ensemble size is impractical and prevents the use of EKI in high-dimensional problems. ‘Localization’ has been used for many years in EnKF to break the subspace property in a way that a localized EnKF can solve high-dimensional problems with a modest ensemble size, independently of the number of unknowns. Here, we study localization of the EKI and demonstrate how a localized EKI (LEKI) can solve high-dimensional inverse problems with a modest ensemble size. Our analysis is mathematically rigorous and applies to the continuous time limit of the EKI. Specifically, we can prove an intended ensemble collapse and convergence guarantees with an ensemble size that is less than the number of unknowns, which sets this work apart from the current state-of-the-art. We illustrate our theory with numerical experiments where some of our mathematical assumptions may only be approximately valid.

Butterfly Transforms for Efficient Representation of Spatially Variant Point Spread Functions in Bayesian Imaging

Bayesian imaging algorithms are becoming increasingly important in, e.g., astronomy, medicine and biology. Given that many of these algorithms compute iterative solutions to high-dimensional inverse problems, the efficiency and accuracy of the instrument response representation are of high importance for the imaging process. For efficiency reasons, point spread functions, which make up a large fraction of the response functions of telescopes and microscopes, are usually assumed to be spatially invariant in a given field of view and can thus be represented by a convolution. For many instruments, this assumption does not hold and degrades the accuracy of the instrument representation. Here, we discuss the application of butterfly transforms, which are linear neural network structures whose sizes scale sub-quadratically with the number of data points. Butterfly transforms are efficient by design, since they are inspired by the structure of the Cooley-Tukey fast Fourier transform. In this work, we combine them in several ways into butterfly networks, compare the different architectures with respect to their performance and identify a representation that is suitable for the efficient representation of a synthetic spatially variant point spread function up to a 1% error. Furthermore, we show its application in a short synthetic example.

Reliable amortized variational inference with physics-based latent distribution correction

Bayesian inference for high-dimensional inverse problems is computationally costly and requires selecting a suitable prior distribution. Amortized variational inference addresses these challenges by pretraining a neural network that acts as a surrogate conditional distribution that approximates the posterior distribution not only for one instance of the observed data but also for the distribution of the data pertaining to a specific inverse problem. When fed previously unseen data, the neural network — in our case, a conditional normalizing flow — provides the posterior samples at virtually no cost. However, the accuracy of amortized variational inference relies on the availability of high-fidelity training data, which seldom exist in geophysical inverse problems because of the earth’s heterogeneous subsurface. In addition, the network is prone to errors if evaluated over data that are not drawn from the training data distribution. As such, we have aimed to increase the resilience of amortized variational inference in the presence of moderate data distribution shifts. We achieve this via a correction to the conditional normalizing flow’s latent distribution that improves the approximation to the posterior distribution for the data at hand. The correction involves relaxing the standard Gaussian assumption on the latent distribution and parameterizing it via a Gaussian distribution with an unknown mean and (diagonal) covariance. These unknowns are then estimated by minimizing the Kullback-Leibler divergence between the corrected and the (physics-based) true posterior distributions. Although generic and applicable to other inverse problems by means of a linearized seismic imaging example, we find that our correction step improves the robustness of amortized variational inference with respect to changes in the number of seismic sources, noise variance, and shifts in the prior distribution. This approach, given noisy seismic data simulated via the linearized Born modeling, provides a seismic image with limited artifacts and an assessment of its uncertainty at approximately the same cost as five reverse time migrations.

Minimizing the distortions in electrophysiological source imaging of cortical oscillatory activity via Spectral Structured Sparse Bayesian Learning.

Oscillatory processes at all spatial scales and on all frequencies underpin brain function. Electrophysiological Source Imaging (ESI) is the data-driven brain imaging modality that provides the inverse solutions to the source processes of the EEG, MEG, or ECoG data. This study aimed to carry out an ESI of the source cross-spectrum while controlling common distortions of the estimates. As with all ESI-related problems under realistic settings, the main obstacle we faced is a severely ill-conditioned and high-dimensional inverse problem. Therefore, we opted for Bayesian inverse solutions that posited a priori probabilities on the source process. Indeed, rigorously specifying both the likelihoods and a priori probabilities of the problem leads to the proper Bayesian inverse problem of cross-spectral matrices. These inverse solutions are our formal definition for cross-spectral ESI (cESI), which requires a priori of the source cross-spectrum to counter the severe ill-condition and high-dimensionality of matrices. However, inverse solutions for this problem were NP-hard to tackle or approximated within iterations with bad-conditioned matrices in the standard ESI setup. We introduce cESI with a joint a priori probability upon the source cross-spectrum to avoid these problems. cESI inverse solutions are low-dimensional ones for the set of random vector instances and not random matrices. We achieved cESI inverse solutions through the variational approximations via our Spectral Structured Sparse Bayesian Learning (ssSBL) algorithm https://github.com/CCC-members/Spectral-Structured-Sparse-Bayesian-Learning. We compared low-density EEG (10-20 system) ssSBL inverse solutions with reference cESIs for two experiments: (a) high-density MEG that were used to simulate EEG and (b) high-density macaque ECoG that were recorded simultaneously with EEG. The ssSBL resulted in two orders of magnitude with less distortion than the state-of-the-art ESI methods. Our cESI toolbox, including the ssSBL method, is available at https://github.com/CCC-members/BC-VARETA_Toolbox.

Weighted Radial Basis Collocation Method for the Nonlinear Inverse Helmholtz Problems

In this paper, a meshfree weighted radial basis collocation method associated with the Newton’s iteration method is introduced to solve the nonlinear inverse Helmholtz problems for identifying the parameter. All the measurement data can be included in the least-squares solution, which can avoid the iteration calculations for comparing the solutions with part of the measurement data in the Galerkin-based methods. Appropriate weights are imposed on the boundary conditions and measurement conditions to balance the errors, which leads to the high accuracy and optimal convergence for solving the inverse problems. Moreover, it is quite easy to extend the solution process of the one-dimensional inverse problem to high-dimensional inverse problem. Nonlinear numerical examples include one-, two- and three-dimensional inverse Helmholtz problems of constant and varying parameter identification in regular and irregular domains and show the high accuracy and exponential convergence of the presented method.

Retrieval of initial condition for Burgers' equation using reduced-order EnKF via POD-based sparse observations

<p style='text-indent:20px;'>The data assimilation (DA) is a popular method to solve uncertainty quantification problem which has attracted more and more attention in disaster assessment and climate change research. However, the DA usually faces the computational problem associated with "the curse of dimensionality'' when the high-dimensional inverse problem is involved. For this, the usual strategy is to combine the dimensionality reduction method with the DA to mitigate the computational demand in practical application. In this paper, we develop a reduced-order model (ROM)-based EnKF used for retrieving the initial condition, where the prediction model and the observation model are all POD-dependent. In addition to the reduced-order model (ROM) being prediction model, the idea of compressed sensing motivates the acquisition of the observation model that allows an efficient combination of the EnKF with QD algorithm to obtain the optimal sparse observation locations. In this way, computational acceleration is gained in retrieving the initial condition. The effectiveness of this algorithm is demonstrated through retrieving the initial condition of an one-dimensional (1-D) Burgers' equation. Experimental results show that a satisfactory retrieval result is obtained using the current ROM-based EnKF with the computational (CPU) time, <inline-formula><tex-math id="M1">\begin{document}$ 42.64\rm{s} $\end{document}</tex-math></inline-formula>, nearly 25 times faster, and the error result <inline-formula><tex-math id="M2">\begin{document}$ 7.47\times 10^{-3} $\end{document}</tex-math></inline-formula>, two order of magnitude <inline-formula><tex-math id="M3">\begin{document}$ O(10^2) $\end{document}</tex-math></inline-formula> more accurate, than the implementation of the traditional full-order model (FOM)-based EnKF where the computational instability is observed at each iteration though in the same condition as set in the ROM-based EnKF, leading to a poor result (the CPU time <inline-formula><tex-math id="M4">\begin{document}$ 1057.97\rm{s} $\end{document}</tex-math></inline-formula>, and the error accuracy <inline-formula><tex-math id="M5">\begin{document}$ 4.06\times 10^{-1}) $\end{document}</tex-math></inline-formula>. The present study extends the application of ROM to efficiently dealing with ill-posed problem as a regularization strategy through replacing the FOM with the ROM such that the low-order truncation of the POD basis leads to the ROM with as few degrees of freedom as possible, which can efficiently inhibit the higher frequency errors, and hence achieve stability of the proposed approach.</p>

Machine learning enables precise holographic characterization of colloidal materials in real time.

Holographic particle characterization uses in-line holographic video microscopy to track and characterize individual colloidal particles dispersed in their native fluid media. Applications range from fundamental research in statistical physics to product development in biopharmaceuticals and medical diagnostic testing. The information encoded in a hologram can be extracted by fitting to a generative model based on the Lorenz-Mie theory of light scattering. Treating hologram analysis as a high-dimensional inverse problem has been exceptionally successful, with conventional optimization algorithms yielding nanometer precision for a typical particle's position and part-per-thousand precision for its size and index of refraction. Machine learning previously has been used to automate holographic particle characterization by detecting features of interest in multi-particle holograms and estimating the particles' positions and properties for subsequent refinement. This study presents an updated end-to-end neural-network solution called CATCH (Characterizing and Tracking Colloids Holographically) whose predictions are fast, precise, and accurate enough for many real-world high-throughput applications and can reliably bootstrap conventional optimization algorithms for the most demanding applications. The ability of CATCH to learn a representation of Lorenz-Mie theory that fits within a diminutive 200 kB hints at the possibility of developing a greatly simplified formulation of light scattering by small objects.

SOLVING HIGH-DIMENSIONAL INVERSE PROBLEMS WITH AUXILIARY UNCERTAINTY VIA OPERATOR LEARNING WITH LIMITED DATA

In complex large-scale systems such as climate, important effects are caused by a combination of confounding processes that are not fully observable. The identification of sources from observations of the system state is vital for attribution and prediction, which inform critical policy decisions. The difficulty of these types of inverse problems lies in the inability to isolate sources and the cost of simulating computational models. Surrogate models may enable the many-query algorithms required for source identification, but data challenges arise from high dimensionality of the state and source, limited ensembles of costly model simulations to train a surrogate model, and few and potentially noisy state observations for inversion due to measurement limitations. The influence of auxiliary processes adds an additional layer of uncertainty that further confounds source identification. We introduce a framework based on (1) calibrating deep neural network surrogates to the flow maps provided by an ensemble of simulations obtained by varying sources, and (2) using these surrogates in a Bayesian framework to identify sources from observations via optimization. Focusing on an atmospheric dispersion exemplar, we find that the expressive and computationally efficient nature of the deep neural network operator surrogates in appropriately reduced dimension allows for source identification with uncertainty quantification using limited data. Introducing a variable wind field as an auxiliary process, we find that a Bayesian approximation error approach is essential for reliable source inversion when uncertainty due to wind stresses the algorithm.

Implicit regularization with strongly convex bias: Stability and acceleration

Implicit regularization refers to the property of optimization algorithms to be biased towards a certain class of solutions. This property is relevant to understand the behavior of modern machine learning algorithms as well as to design efficient computational methods. While the case where the bias is given by a Euclidean norm is well understood, implicit regularization schemes for more general classes of biases are much less studied. In this work, we consider the case where the bias is given by a strongly convex functional, in the context of linear models, and data possibly corrupted by noise. In particular, we propose and analyze accelerated optimization methods and highlight a trade-off between convergence speed and stability. Theoretical findings are complemented by an empirical analysis on high-dimensional inverse problems in machine learning and signal processing, showing excellent results compared to the state of the art.

Proximal nested sampling for high-dimensional Bayesian model selection

Bayesian model selection provides a powerful framework for objectively comparing models directly from observed data, without reference to ground truth data. However, Bayesian model selection requires the computation of the marginal likelihood (model evidence), which is computationally challenging, prohibiting its use in many high-dimensional Bayesian inverse problems. With Bayesian imaging applications in mind, in this work we present the proximal nested sampling methodology to objectively compare alternative Bayesian imaging models for applications that use images to inform decisions under uncertainty. The methodology is based on nested sampling, a Monte Carlo approach specialised for model comparison, and exploits proximal Markov chain Monte Carlo techniques to scale efficiently to large problems and to tackle models that are log-concave and not necessarily smooth (e.g., involving ell _1 or total-variation priors). The proposed approach can be applied computationally to problems of dimension {mathcal {O}}(10^6) and beyond, making it suitable for high-dimensional inverse imaging problems. It is validated on large Gaussian models, for which the likelihood is available analytically, and subsequently illustrated on a range of imaging problems where it is used to analyse different choices of dictionary and measurement model.

Fast inverse estimation of hydraulic conductivity field based on a deep convolutional-cycle generative adversarial neural network

This study proposes a novel deep convolutional-cycle generative adversarial neural network (DC-CGAN) for estimating the hydraulic conductivity field K from the observed head field H. The estimation of the heterogeneous hydraulic conductivity is apparently a high-dimensional inverse problem. Although the indirect approach is popularly utilized to solve the inverse problem, when compared with the indirect approach, the direct approach can conduct the inversion task faster without any iteration process. The key of the direct approach is to establish the inverse operator which maps the observed head to hydraulic conductivity. However, the complex inverse mapping relationship and high-dimensional estimation have always challenged the direct approach. Owing to the potential ability to handle high-dimensional images, six residual blocks with convolutional kernels have been designed as the generators of DC-CGAN to capture the potential distribution for enhancing the generation of estimation images similar to the real hydraulic conductivity. However, the traditional adversarial pattern of one-way generator and discriminator might suffer from the equifinality from different parameters. Therefore, to alleviate the equifinality, a cycle adversarial training pattern involving two generators and discriminators has been employed to supervise the map. Particularly, this pattern performs a forward map from H to K and then the reconstructed K back to H. The reconstruction process has been a supervision to the forward map process. A hypothetical case with 40×30 unknown K values has been designed to evaluate the performance of DC-CGAN, and the results indicate that DC-CGAN has achieved a desired generalization accuracy of MRE of 9.25 % and SSIM of 0.85 with the swift running speed of 0.2 s, which outperforms indirect approach. The application of the proposed DC-CGAN shows that it can provide a reliable and fast high-dimensional estimation task for the heterogenous hydraulic conductivity.

Hydrogeological multiple-point statistics inversion by adaptive sequential Monte Carlo

For strongly non-linear and high-dimensional inverse problems, Markov chain Monte Carlo (MCMC) methods may fail to properly explore the posterior probability density function (PDF) given a realistic computational budget and are generally poorly amenable to parallelization. Particle methods approximate the posterior PDF using the states and weights of a population of evolving particles and they are very well suited to parallelization. We focus on adaptive sequential Monte Carlo (ASMC), an extension of annealed importance sampling (AIS). In AIS and ASMC, importance sampling is performed over a sequence of intermediate distributions, known as power posteriors, linking the prior to the posterior PDF. The AIS and ASMC algorithms also provide estimates of the evidence (marginal likelihood) as needed for Bayesian model selection, at basically no additional cost. ASMC performs better than AIS as it adaptively tunes the tempering schedule and performs resampling of particles when the variance of the particle weights becomes too large. We consider a challenging synthetic groundwater transport inverse problem with a categorical channelized 2D hydraulic conductivity field defined such that the posterior facies distribution includes two distinct modes. The model proposals are obtained by iteratively re-simulating a fraction of the current model using conditional multiple-point statistics (MPS) simulations. We examine how ASMC explores the posterior PDF and compare with results obtained with parallel tempering (PT), a state-of-the-art MCMC inversion approach that runs multiple interacting chains targeting different power posteriors. For a similar computational budget, ASMC outperforms PT as the ASMC-derived models fit the data better and recover the reference likelihood. Moreover, we show that ASMC partly retrieves both posterior modes, while none of them is recovered by PT. Lastly, we demonstrate how the power posteriors obtained by ASMC can be used to assess the influence of the assumed data errors on the posterior means and variances, as well as on the evidence. We suggest that ASMC can advantageously replace MCMC for solving many challenging inverse problems arising in the field of water resources.

α-deep Probabilistic Inference (α-DPI): Efficient Uncertainty Quantification from Exoplanet Astrometry to Black Hole Feature Extraction

Inference is crucial in modern astronomical research, where hidden astrophysical features and patterns are often estimated from indirect and noisy measurements. Inferring the posterior of hidden features, conditioned on the observed measurements, is essential for understanding the uncertainty of results and downstream scientific interpretations. Traditional approaches for posterior estimation include sampling-based methods and variational inference (VI). However, sampling-based methods are typically slow for high-dimensional inverse problems, while VI often lacks estimation accuracy. In this paper, we propose α-deep probabilistic inference, a deep learning framework that first learns an approximate posterior using α-divergence VI paired with a generative neural network, and then produces more accurate posterior samples through importance reweighting of the network samples. It inherits strengths from both sampling and VI methods: it is fast, accurate, and more scalable to high-dimensional problems than conventional sampling-based approaches. We apply our approach to two high-impact astronomical inference problems using real data: exoplanet astrometry and black hole feature extraction.