Abstract: United States. Dept. of Energy. Office of Advanced Scientific Computing Research (Grant DE-SC0003908)

Abstract: We propose and analyze a new formulation of the Linear Sampling Method that uses an exact characterization of the targets shape in terms of the so-called farfield operator (at a fixed frequency). This characterization is based on constructing nearby solutions of the farfield equation using minimizing sequences of a least squares cost functional with an appropriate penalty term. We first provide a general framework for the theoretical foundation of the method in the case of noise-free and noisy measurements operator. We then explicit applications for the case of inhomogeneous inclusions and indicate possible straightforward generalizations. We finally validate the method through some numerical tests and compare the performances with classical LSM and the factorization methods.

Abstract: In this paper, we consider the inverse problem of determining the permeability of the subsurface from hydraulic head measurements, within the framework of a steady Darcy model of groundwater flow. We study geometrically defined prior permeability fields, which admit layered, fault and channel structures, in order to mimic realistic subsurface features; within each layer we adopt either a constant or continuous function representation of the permeability. This prior model leads to a parameter identification problem for a finite number of unknown parameters determining the geometry, together with either a finite number of permeability values (in the constant case) or a finite number of fields (in the continuous function case). We adopt a Bayesian framework showing the existence and well-posedness of the posterior distribution. We also introduce novel Markov chain Monte Carlo (MCMC) methods, which exploit the different character of the geometric and permeability parameters, and build on recent advances in function space MCMC. These algorithms provide rigorous estimates of the permeability, as well as the uncertainty associated with it, and only require forward model evaluations. No adjoint solvers are required and hence the methodology is applicable to black-box forward models. We then use these methods to explore the posterior and to illustrate the methodology with numerical experiments.

Abstract: Airborne time-domain electromagnetic surveys are effective tools for mineral exploration and geologic mapping. 3D inversion of airborne electromagnetic data is a challenging computational problem. The size of the surveys and the spatial resolution required to adequately discretize the transmitters and receivers results in very large meshes. Solving the forward problem repeatedly on such a mesh can quickly become impractical. Fortunately, using a single mesh for both the forward and inverse problem for all of the transmitters is not necessary. The forward problem for a single source or a small group of sources can be solved on different meshes, each of which only needs to be locally refined close to the selected transmitters and receivers. Away from the selected transmitters and receivers, the mesh can be coarsened. The forward problem can then be broken into a number of highly parallel problems. Each forward modelling mesh is optimized specifically to the selected transmitters and receivers and has far fewer cells than the fine inversion mesh. Further efficiency can be gained by using stochastic Gauss–Newton methods where a stochastic approximation to the gradient, Hessian or both are used. In this paper, we present new algorithms for airborne data inversion and their implementation using a finite volume discretization on OcTree meshes. We demonstrate our approach on a large-scale synthetic versatile time-domain electromagnetic surveying data set.

Abstract: Solving inverse problems based on computationally demanding forward models is ubiquitously difficult since one is necessarily limited to just a few observations of the response surface. The usual practice is to replace the response surface with a surrogate. However, this approach induces additional uncertainties on the posterior distributions. The main contribution of this work is the reformulation of the Bayesian solution of the inverse problem when the expensive forward model is replaced by the surrogate. We derive three approximations of the reformulated solution with increasing complexity and fidelity. We demonstrate numerically that the proposed approximations capture theuncertaintyofthesolutionoftheinverseprobleminducedbythefactthatthe forward model is replaced by a finite number of simulations. We demonstrate our approach in two different problems: locating the contamination source of a diffusive process and inferring the permeability field of an oil reservoir based on measurements of the oil-cut curves.

Abstract: We consider the Riemann manifold Hamiltonian Monte Carlo (RMHMC) method for solving statistical inverse problems governed by partial differential equations (PDEs). The Bayesian framework is employed to cast the inverse problem into the task of statistical inference whose solution is the posterior distribution in infinite dimensional parameter space conditional upon observation data and Gaussian prior measure. We discretize both the likelihood and the prior using the H1-conforming finite element method together with a matrix transfer technique. The power of the RMHMC method is that it exploits the geometric structure induced by the PDE constraints of the underlying inverse problem. Consequently, each RMHMC posterior sample is almost uncorrelated/independent from the others providing statistically efficient Markov chain simulation. However this statistical efficiency comes at a computational cost. This motivates us to consider computationally more efficient strategies for RMHMC. At the heart of our construction is the fact that for Gaussian error structures the Fisher information matrix coincides with the Gauss–Newton Hessian. We exploit this fact in considering a computationally simplified RMHMC method combining state-of-the-art adjoint techniques and the superiority of the RMHMC method. Specifically, we first form the Gauss–Newton Hessian at the maximum a posteriori point and then use it as a fixed constant metric tensor throughout RMHMC simulation. This eliminates the need for the computationally costly differential geometric Christoffel symbols, which in turn greatly reduces computational effort at a corresponding loss of sampling efficiency. We further reduce the cost of forming the Fisher information matrix by using a low rank approximation via a randomized singular value decomposition technique. This is efficient since a small number of Hessian-vector products are required. The Hessian-vector product in turn requires only two extra PDE solves using the adjoint technique. Various numerical results up to 1025 parameters are presented to demonstrate the ability of the RMHMC method in exploring the geometric structure of the problem to propose (almost) uncorrelated/independent samples that are far away from each other, and yet the acceptance rate is almost unity. The results also suggest that for the PDE models considered the proposed fixed metric RMHMC can attain almost as high a quality performance as the original RMHMC, i.e. generating (almost) uncorrelated/independent samples, while being two orders of magnitude less computationally expensive.

Abstract: The aim of this paper is to test and analyse a novel technique for image reconstruction in positron emission tomography, which is based on (total variation) regularization on both the image space and the projection space. We formulate our variational problem considering both total variation penalty terms on the image and on an idealized sinogram to be reconstructed from a given Poisson distributed noisy sinogram. We prove existence, uniqueness and stability results for the proposed model and provide some analytical insight into the structures favoured by joint regularization. For the numerical solution of the corresponding discretized problem we employ the split Bregman algorithm and extensively test the approach in comparison to standard total variation regularization on the image. The numerical results show that an additional penalty on the sinogram performs better on reconstructing images with thin structures.

Abstract: This work investigates the electrical impedance tomography (EIT) problem in the case when only one or two pairs of Cauchy data is available, which is known to be very difficult in achieving high reconstruction quality owing to its severely ill-posed nature. We propose a simple and efficient direct sampling method (DSM) to locate inhomogeneous inclusions. A new probing function based on the dipole potential is introduced to construct an indicator function for imaging the inclusions. Explicit formulae for the probing and indicator functions are derived in the case when the sampling domain is of spherical geometry in (n = 2, 3). This new method is easy to implement and computationally cheap. Numerical experiments are presented to demonstrate the robustness and effectiveness of the DSM, which provides a new numerical approach for solving the EIT problem.

Abstract: A frequent matter of debate in Bayesian inversion is the question of which of the two principal point-estimators, the maximum a posteriori (MAP) or the conditional mean (CM) estimate, is to be preferred. As the MAP estimate corresponds to the solution given by variational regularization techniques, this is also a constant matter of debate between the two research areas. Following a theoretical argument—the Bayes cost formalism—the CM estimate is classically preferred for being the Bayes estimator for the mean squared error cost, while the MAP estimate is classically discredited for being only asymptotically the Bayes estimator for the uniform cost function. In this article we present recent theoretical and computational observations that challenge this point of view, in particular for high-dimensional sparsity-promoting Bayesian inversion. Using Bregman distances, we present new, proper convex Bayes cost functions for which the MAP estimator is the Bayes estimator. We complement this finding with results that correct further common misconceptions about MAP estimates. In total, we aim to rehabilitate MAP estimates in linear inverse problems with log-concave priors as proper Bayes estimators.

Abstract: We study a new class of Radon transforms defined on circular cones called the conical Radon transform. In it maps a function to its surface integrals over circular cones, and in it maps a function to its integrals along two rays with a common vertex. Such transforms appear in various mathematical models arising in medical imaging, nuclear industry and homeland security. This paper contains new results about inversion of conical Radon transform with fixed opening angle and vertical central axis in and . New simple explicit inversion formulae are presented in these cases. Numerical simulations were performed to demonstrate the efficiency of the suggested algorithm in 2D.

Abstract: Sparsity optimization for linear least squares problems formulated as non-smooth regularization problems are considered in infinite-dimensional sequence spaces lp with p ∈ [0, 1]. Necessary optimality conditions in the format of a complementarity system are obtained. A monotonically convergent scheme is developed for the case p ∈ (0, 1]. For the case p = 0 a primal dual active set strategy based on the Lagrange multiplier rule is proposed and analyzed for special cases.

Abstract: Inversion of Radon transforms is the mathematical foundation of many modern tomographic imaging modalities. In this paper, we study a conical Radon transform, which is important for computed tomography taking Compton scattering into account. This transform integrates a function in over all conical surfaces having vertices on a hyperplane and symmetry axis orthogonal to this plane. As the main result, we derive exact reconstruction formulas of the filtered back-projection type for inverting this transform. Further, a Fourier slice identity is established.

Abstract: We establish posterior sparsity in Bayesian inversion for systems governed by operator equations with distributed parameter uncertainty subject to noisy observation data δ. We generalize the results and algorithms introduced in C Schillings and C Schwab (2013 Inverse Problems 29 065011) for the particular case of scalar diffusion problems with random coefficients to broad classes of forward problems, including general elliptic and parabolic operators with uncertain coefficients, and in random domains. For countably parametric, deterministic representations of uncertain parameters in the forward problem, which belong to a specified sparsity class, we quantify analytic regularity of the likewise countably parametric, deterministic Bayesian posterior density with respect to a uniform prior on the uncertain parameter sequences and prove that the parametric, deterministic density of the Bayesian posterior belongs to the same sparsity class. Generalizing C Schillings and C Schwab (2013 Inverse Problems 29 065011) and C Schwab and A M Stuart (2012 Inverse Problems 28 045003) the forward problems are converted to countably parametric, deterministic operator equations. Computational Bayesian inversion amounts to numerically evaluating expectations of quantities of interest (QoIs) under the Bayesian posterior, conditional on noisy observation data. Our results imply, on the one hand, sparsity of Legendre (generalized) polynomial chaos expansions of the density of the Bayesian posterior with respect to uniform prior and, on the other hand, convergence rates for data-adaptive Smolyak integration algorithms for computational Bayesian estimation, which are independent of the dimension of the parameter space. We prove, mathematically and computationally, that for uncertain inputs with sufficient sparsity convergence rates are, in particular, superior to Markov chain Monte-Carlo sampling of the posterior, in terms of the number N of instances of the parametric forward problem to be solved.

Abstract: We consider the inverse problem of the reconstruction of the spatially distributed dielectric constant epsilon(r)(x), x is an element of R-3, which is an unknown coefficient in the Maxwell's equations, from time-dependent backscattering experimental radar data associated with a single source of electric pulses. The refractive index is n(x) = root epsilon(r)(x). The coefficient epsilon(r)(x) is reconstructed using a two-stage reconstruction procedure. In the first stage an approximately globally convergent method proposed is applied to get a good first approximation of the exact solution. In the second stage a locally convergent adaptive finite element method is applied, taking the solution of the first stage as the starting point of the minimization of the Tikhonov functional. This functional is minimized on a sequence of locally refined meshes. It is shown here that all three components of interest of targets can be simultaneously accurately imaged: refractive indices, shapes and locations.

Abstract: An inverse problem is called dynamic if the object changes during the data acquisition process. This occurs e.g. in medical applications when fast moving organs like the lungs or the heart are imaged. Most regularization methods are based on the assumption that the object is static during the measuring procedure. Hence, their application in the dynamic case often leads to serious motion artefacts in the reconstruction. Therefore, an algorithm has to take into account the temporal changes of the investigated object. In this paper, a reconstruction method that compensates for the motion of the object is derived for dynamic linear inverse problems. The algorithm is validated at numerical examples from computerized tomography.

Abstract: The topological derivative-based non-iterative imaging algorithm has demonstrated its applicability in limited-aperture inverse scattering problems. However, this has been confirmed through many experimental simulation results, and the reason behind this applicability has not been satisfactorily explained. In this paper, we identify the mathematical structure and certain properties of topological derivatives for the imaging of two-dimensional crack-like thin penetrable electromagnetic inhomogeneities that are completely embedded in a homogeneous material. To this end, we establish a relationship with an infinite series of Bessel functions of integer order of the first kind. Based on the derived structure, we discover a necessary condition for applying topological derivatives in limited-aperture inverse scattering problems, and thus confirm why topological derivatives can be applied. Furthermore, we analyze the structure of multi-frequency topological derivative, and identify why this improves the single-frequency topological derivative in limited-aperture inverse scattering problems. Various numerical simulations are conducted with noisy data, and the results support the derived structure and exhibit certain properties of single- and multi-frequency topological derivatives.

Abstract: We consider the problem of phase retrieval in the Fresnel regime. In recent years, several techniques have been used to solve this problem applying different a priori assumptions on the two-dimensional object in space, such as positivity, finite support and amplitude constraints. In this paper, we propose a new constraint, namely the assumption that the object possesses a sparse representation in a shearlet frame. We show how a shearlet soft-thresholding procedure can be used for phase reconstruction with Fresnel data. As it turns out, the shearlet sparsity constraint yields reconstruction results that are far superior to the support constraint, and comparable to the support plus positivity constraint.

Abstract: In this paper we consider the well-established problem of time differences of arrival (TDOA)-based source localization and propose a comprehensive analysis of its solution for arbitrary sensor measurement and placement. More specifically, we define the TDOA map from the physical space of source locations to the space of range measurements (TDOAs), in the specific case of three receivers in 2D space. We then study the identifiability of the model, giving a complete analytical characterization of the image of this map and its invertibility. This analysis has been conducted in a completely mathematical fashion, using many different tools which make it valid for every sensor configuration. These results are the first step toward the solution of more general problems involving, for example, a larger number of sensors, uncertainty in their placement, or lack of synchronization.

Abstract: We present a finite element analysis of electrical impedance tomography for reconstructing the conductivity distribution from electrode voltage measurements by means of Tikhonov regularization. Two popular choices of the penalty term, i.e., the H1(Ω)-norm smoothness penalty and total variation seminorm penalty, are considered. A piecewise linear finite element method is employed for discretizing the forward model, i.e., the complete electrode model, the conductivity, and the penalty functional. The convergence of the finite element approximations for the Tikhonov model on both polyhedral and smooth curved domains is established. This provides rigorous justifications for the ad hoc discretization procedures. Numerical experiments confirm the convergence analysis.

Abstract: In this paper, we consider an inverse source problem for elliptic partial differential equations with Dirichlet and Neumann boundary data. The unknown source term is to be determined from additional boundary conditions. Unlike the existing methods found in the literature, which usually use some of the boundary conditions to form a boundary value problem for the elliptic partial differential equation and the remaining boundary conditions in the objective functional for optimization to determine the source term, the novel method that we propose here has coupled complex boundary conditions. We use a complex elliptic partial differential equation with a Robin boundary condition coupling the Dirichlet and Neumann boundary data, and optimize with respect to the imaginary part of the solution in the domain to determine the source term. Then, on the basis of the complex boundary value problem, Tikhonov regularization is used to obtain a stable approximate source function and the finite element method is used for discretization. Theoretical analysis is given for both the continuous model and the discrete model. Several numerical examples are provided to show the usefulness of the proposed coupled complex boundary method.

Abstract: Variable splitting schemes for the function space version of the image reconstruction problem with total variation regularization (TV-problem) in its primal and pre-dual formulations are considered. For the primal splitting formulation, while existence of a solution cannot be guaranteed, it is shown that quasi-minimizers of the penalized problem are asymptotically related to the solution of the original TV-problem. On the other hand, for the pre-dual formulation, a family of parametrized problems is introduced and a parameter dependent contraction of an associated fixed point iteration is established. Moreover, the theory is validated by numerical tests. Additionally, the augmented Lagrangian approach is studied, details on an implementation on a staggered grid are provided and numerical tests are shown.

Abstract: The history of Bayesian methods dates back to the original works of Reverend Thomas Bayes and Pierre-Simon Laplace: the former laid down some of the basic principles on inverse probability in his classic article ‘An essay towards solving a problem in the doctrine of chances’ that was read posthumously in the Royal Society in 1763. Laplace, on the other hand, in his ‘Memoirs on inverse probability’ of 1774 developed the idea of updating beliefs and wrote down the celebrated Bayes’ formula in the form we know today. Although not identified yet as a framework for investigating inverse problems, Laplace used the formalism very much in the spirit it is used today in the context of inverse problems, e.g., in his study of the distribution of comets. With the evolution of computational tools, Bayesian methods have become increasingly popular in all fields of human knowledge in which conclusions need to be drawn based on incomplete and noisy data. Needless to say, inverse problems, almost by definition, fall into this category. Systematic work for developing a Bayesian inverse problem framework can arguably be traced back to the 1980s, see [25] and [26] (the original first edition being published by Elsevier in 1987), although articles on Bayesian methodology applied to inverse problems, in particular in geophysics, had appeared much earlier. For more recent overviews, see [4, 15, 24]. Today, as testified by the articles in this special issue, the Bayesian methodology as a framework for considering inverse problems has gained a lot of popularity, and it has integrated very successfully with many traditional inverse problems ideas and techniques, providing novel ways to interpret and implement traditional procedures in numerical analysis, computational statistics, signal analysis and data assimilation. The range of applications where the Bayesian framework has been fundamental goes from geophysics, engineering and imaging to astronomy, life sciences and economy, and continues to grow. There is no question that Bayesian methodology constitutes one of the cornerstones in the increasing research area of uncertainty quantification, a concept that is right in the heart of Bayesian theory. This special issue widens the spectrum of Bayesian inverse problems to include, in addition to Bayesian methodological developments designed for a variety of applications known already to the readership of the journal Inverse Problems, fields that have been less represented in this journal and scientific conferences around inverse problems. In the Bayesian framework, the complete solution of the inverse problem is the posterior distribution of the unknown, and in order to explore the posterior distribution, effective sampling strategies need to be developed. The standard versions of Markov Chain Monte Carlo (MCMC) methods such as the Metropolis–Hastings algorithm are rather easy to describe, but efficient numerical implementation is often a problem for various reasons. The article [2] discusses the implementation strategies of a particular version of MCMC, the Riemannian manifold Hamiltonian Monte Carlo algorithm, with an application to the steady state inverse heat conduction problem. The article [9] addresses a trans-dimensional problem, or a problem in which the dimension of the inverse problem is not known, requiring a Inverse Problems

Abstract: We prove a global Holder stability estimate for a hybrid inverse problem combining microwave imaging and ultrasound. The principal features of this result are that we assume to have access to measurements associated with a single, arbitrary and possibly sign changing solution of a Schrodinger equation, and that zero is allowed to be an eigenvalue of the equation.

Abstract: This paper concerns the reconstruction of a fully anisotropic conductivity tensor γ from internal current densities of the form J = γ∇u, where u solves a second-order elliptic equation ∇ (γ∇u) = 0 on a bounded domain X with prescribed boundary conditions. A minimum number of n + 2 such functionals known on Y⊂X, where n is the spatial dimension, is sufficient to guarantee a unique and explicit reconstruction of γ locally on Y. Moreover, we show that γ is reconstructed with a loss of one derivative compared to errors in the measurement of J in the general case and no loss of derivatives in the special case where γ is scalar. We also describe linear combinations of mixed partial derivatives of γ that exhibit better stability properties and hence can be reconstructed with better resolution in practice.

Abstract: We propose a novel numerical inversion algorithm for the coefficients of parabolic partial differential equations, based on model reduction. The study is motivated by the application of controlled source electromagnetic exploration, where the unknown is the subsurface electrical resistivity and the data are time resolved surface measurements of the magneticfield. The algorithm presented in this paper considers inversion in one and two dimensions. The reduced model is obtained with rational interpolation in the frequency (Laplace) domain and a rational Krylov subspace projection method. It amounts to a nonlinear mapping from the function space of the unknown resistivity to the small dimensional space of the parameters of the reduced model. We use this mapping as a nonlinear preconditioner for the Gauss‐Newton iterative solution of the inverse problem. The advantage of the inversion algorithm is twofold. First, the nonlinear preconditioner resolves most of the nonlinearity of the problem. Thus the iterations are less likely to get stuck in local minima and the convergence is fast. Second, the inversion is computationally efficient because it avoids repeated accurate simulations of the time-domain response. We study the stability of the inversion algorithm for various rational Krylov subspaces, and assess its performance with numerical experiments.

Abstract: In this paper, we develop a novel Bayesian approach to the problem of estimating neural currents in the brain from a fixed distribution of magnetic field (called topography), measured by magnetoencephalography. Differently from recent studies that describe inversion techniques, such as spatio-temporal regularization/filtering, in which neural dynamics always plays a role, we face here a purely static inverse problem. Neural currents are modelled as an unknown number of current dipoles, whose state space is described in terms of a variable-dimension model. Within the resulting Bayesian framework, we set up a sequential Monte Carlo sampler to explore the posterior distribution. An adaptation technique is employed in order to effectively balance the computational cost and the quality of the sample approximation. Then, both the number and the parameters of the unknown current dipoles are simultaneously estimated. The performance of the method is assessed by means of synthetic data, generated by source configurations containing up to four dipoles. Eventually, we describe the results obtained by analysing data from a real experiment, involving somatosensory evoked fields, and compare them to those provided by three other methods.

Abstract: In this paper we establish a factorization method for recovering the location and shape of an acoustic bounded obstacle with using the near-field data, corresponding to infinitely many incident point sources. The obstacle is allowed to be an impenetrable scatterer of sound-soft, sound-hard or impedance type or a penetrable scatterer. An outgoing-to-incoming operator is constructed for facilitating the factorization of the near-field operator, which can be easily implemented numerically. Numerical examples are presented to demonstrate the feasibility and effectiveness of our inversion algorithm, including the case where limited aperture near-field data are available only.

Abstract: Within the framework of linear elasticity we assume the availability of internal full-field measurements of the continuum deformations of a non-homogeneous isotropic solid. The aim is the quantitative reconstruction of the associated moduli. A simple gradient system for the sought constitutive parameters is derived algebraically from the momentum equation, whose coefficients are expressed in terms of the measured displacement fields and their spatial derivatives. Direct integration of this system is discussed to finally demonstrate the inexpediency of such an approach when dealing with noisy data. Upon using polluted measurements, an alternative variational formulation is deployed to invert for the physical parameters. Analysis of this latter inversion procedure provides existence and uniqueness results while the reconstruction stability with respect to the measurements is investigated. As the inversion procedure requires differentiating the measurements twice, a numerical differentiation scheme based on an ad hoc regularization then allows an optimally stable reconstruction of the sought moduli. Numerical results are included to illustrate and assess the performance of the overall approach.

Abstract: We consider the inverse problem of determining the Lame moduli for a piecewise constant elasticity tensor , where is a known finite partition of the body Ω, from the Dirichlet-to-Neumann map. We prove that Lipschitz stability estimates can be derived under regularity assumptions on the interfaces.