Search

We discuss solution algorithms for calibrating a tumor growth model using imaging data posed as a deterministic inverse problem. The forward model consists of a nonlinear and time-dependent reaction-diffusion partial differential equation (PDE) with unknown parameters (diffusivity and proliferation rate) being spatial fields. We use a dimension-independent globalized, inexact Newton Conjugate Gradient algorithm to solve the PDE-constrained optimization. The required gradient and Hessian actions are also presented using the adjoint method and Lagrangian formalism.

Ultrasound computed tomography (USCT) is an emerging medical imaging modality that holds great promise for improving human health. Full-waveform inversion (FWI)-based image reconstruction methods account for the relevant wave physics to produce high spatial resolution images of the acoustic properties of the breast tissues. A practical USCT design employs a circular ring-array comprised of elevation-focused ultrasonic transducers, and volumentric imaging is achieved by translating the ring-array orthogonally to the imaging plane. In commonly deployed slice-by-slice (SBS) reconstruction approaches, the three-dimensional (3D) volume is reconstructed by stacking together two-dimensional (2D) images reconstructed for each position of the ring-array. A limitation of the SBS reconstruction approach is that it does not account for 3D wave propagation physics and the focusing properties of the transducers, which can result in significant image artifacts and inaccuracies. To perform 3D image reconstruction when elevation-focused transducers are employed, a numerical description of the focusing properties of the transducers should be included in the forward model. To address this, a 3D computational model of an elevation-focused transducer is developed to enable 3D FWI-based reconstruction methods to be deployed in ring-array-based USCT. The focusing is achieved by applying a spatially varying temporal delay to the ultrasound pulse (emitter mode) and recorded signal (receiver mode). The proposed numerical transducer model is quantitatively validated and employed in computer-simulation studies that demonstrate its use in image reconstruction for ring-array USCT

Medical imaging systems are often evaluated and optimized via objective, or task-specific, measures of image quality (IQ) that quantify the performance of an observer on a specific clinically-relevant task. The performance of the Bayesian Ideal Observer (IO) sets an upper limit among all observers, numerical or human, and has been advocated for use as a figure-of-merit (FOM) for evaluating and optimizing medical imaging systems. However, the IO test statistic corresponds to the likelihood ratio that is intractable to compute in the majority of cases. A sampling-based method that employs Markov-Chain Monte Carlo (MCMC) techniques was previously proposed to estimate the IO performance. However, current applications of MCMC methods for IO approximation have been limited to a small number of situations where the considered distribution of to-be-imaged objects can be described by a relatively simple stochastic object model (SOM). As such, there remains an important need to extend the domain of applicability of MCMC methods to address a large variety of scenarios where IO-based assessments are needed but the associated SOMs have not been available. In this study, a novel MCMC method that employs a generative adversarial network (GAN)-based SOM, referred to as MCMC-GAN, is described and evaluated. The MCMC-GAN method was quantitatively validated by use of test-cases for which reference solutions were available. The results demonstrate that the MCMC-GAN method can extend the domain of applicability of MCMC methods for conducting IO analyses of medical imaging systems., Comment: Submitted to IEEE Transactions on Medical Imaging

Dynamic imaging is essential for analyzing various biological systems and behaviors but faces two main challenges: data incompleteness and computational burden. For many imaging systems, high frame rates and short acquisition times require severe undersampling, which leads to data incompleteness. Multiple images may then be compatible with the data, thus requiring special techniques (regularization) to ensure the uniqueness of the reconstruction. Computational and memory requirements are particularly burdensome for three-dimensional dynamic imaging applications requiring high resolution in both space and time. Exploiting redundancies in the object's spatiotemporal features is key to addressing both challenges. This contribution investigates neural fields, or implicit neural representations, to model the sought-after dynamic object. Neural fields are a particular class of neural networks that represent the dynamic object as a continuous function of space and time, thus avoiding the burden of storing a full resolution image at each time frame. Neural field representation thus reduces the image reconstruction problem to estimating the network parameters via a nonlinear optimization problem (training). Once trained, the neural field can be evaluated at arbitrary locations in space and time, allowing for high-resolution rendering of the object. Key advantages of the proposed approach are that neural fields automatically learn and exploit redundancies in the sought-after object to both regularize the reconstruction and significantly reduce memory storage requirements. The feasibility of the proposed framework is illustrated with an application to dynamic image reconstruction from severely undersampled circular Radon transform data., 13 pages, 10 figures

Tomographic imaging is in general an ill-posed inverse problem. Typically, a single regularized image estimate of the sought-after object is obtained from tomographic measurements. However, there may be multiple objects that are all consistent with the same measurement data. The ability to generate such alternate solutions is important because it may enable new assessments of imaging systems. In principle, this can be achieved by means of posterior sampling methods. In recent years, deep neural networks have been employed for posterior sampling with promising results. However, such methods are not yet for use with large-scale tomographic imaging applications. On the other hand, empirical sampling methods may be computationally feasible for large-scale imaging systems and enable uncertainty quantification for practical applications. Empirical sampling involves solving a regularized inverse problem within a stochastic optimization framework to obtain alternate data-consistent solutions. In this work, a new empirical sampling method is proposed that computes multiple solutions of a tomographic inverse problem that are consistent with the same acquired measurement data. The method operates by repeatedly solving an optimization problem in the latent space of a style-based generative adversarial network (StyleGAN), and was inspired by the Photo Upsampling via Latent Space Exploration (PULSE) method that was developed for super-resolution tasks. The proposed method is demonstrated and analyzed via numerical studies that involve two stylized tomographic imaging modalities. These studies establish the ability of the method to perform efficient empirical sampling and uncertainty quantification., Comment: Submitted to IEEE Transactions on Medical Imaging

We propose derivative-informed neural operators (DINOs), a general family of neural networks to approximate operators as infinite-dimensional mappings from input function spaces to output function spaces or quantities of interest. After discretizations both inputs and outputs are high-dimensional. We aim to approximate not only the operators with improved accuracy but also their derivatives (Jacobians) with respect to the input function-valued parameter to empower derivative-based algorithms in many applications, e.g., Bayesian inverse problems, optimization under parameter uncertainty, and optimal experimental design. The major difficulties include the computational cost of generating derivative training data and the high dimensionality of the problem leading to large training cost. To address these challenges, we exploit the intrinsic low-dimensionality of the derivatives and develop algorithms for compressing derivative information and efficiently imposing it in neural operator training yielding derivative-informed neural operators. We demonstrate that these advances can significantly reduce the costs of both data generation and training for large classes of problems (e.g., nonlinear steady state parametric PDE maps), making the costs marginal or comparable to the costs without using derivatives, and in particular independent of the discretization dimension of the input and output functions. Moreover, we show that the proposed DINO achieves significantly higher accuracy than neural operators trained without derivative information, for both function approximation and derivative approximation (e.g., Gauss-Newton Hessian), especially when the training data are limited.

The Alternating Direction Method of Multipliers (ADMM) provides a natural way of solving inverse problems with multiple partial differential equations (PDE) forward models and nonsmooth regularization. ADMM allows splitting these large-scale inverse problems into smaller, simpler sub-problems, for which computationally efficient solvers are available. In particular, we apply large-scale second-order optimization methods to solve the fully-decoupled Tikhonov regularized inverse problems stemming from each PDE forward model. We use fast proximal methods to handle the nonsmooth regularization term. In this work, we discuss several adaptations (such as the choice of the consensus norm) needed to maintain consistency with the underlining infinite-dimensional problem. We present two imaging applications inspired by electrical impedance tomography and quantitative photoacoustic tomography to demonstrate the proposed method's effectiveness.

Uncertainty quantification of groundwater (GW) aquifer parameters is critical for efficient management and sustainable extraction of GW resources. These uncertainties are introduced by the data, model, and prior information on the parameters. Here we develop a Bayesian inversion framework that uses Interferometric Synthetic Aperture Radar (InSAR) surface deformation data to infer the laterally heterogeneous permeability of a transient linear poroelastic model of a confined GW aquifer. The Bayesian solution of this inverse problem takes the form of a posterior probability density of the permeability. Exploring this posterior using classical Markov chain Monte Carlo (MCMC) methods is computationally prohibitive due to the large dimension of the discretized permeability field and the expense of solving the poroelastic forward problem. However, in many partial differential equation (PDE)-based Bayesian inversion problems, the data are only informative in a few directions in parameter space. For the poroelasticity problem, we prove this property theoretically for a one-dimensional problem and demonstrate it numerically for a three-dimensional aquifer model. We design a generalized preconditioned Crank--Nicolson (gpCN) MCMC method that exploits this intrinsic low dimensionality by using a low-rank based Laplace approximation of the posterior as a proposal, which we build scalably. The feasibility of our approach is demonstrated through a real GW aquifer test in Nevada. The inherently two dimensional nature of InSAR surface deformation data informs a sufficient number of modes of the permeability field to allow detection of major structures within the aquifer, significantly reducing the uncertainty in the pressure and the displacement quantities of interest.

This paper presents the use of element-based algebraic multigrid (AMGe) hierarchies, implemented in the ParELAG (Parallel Element Agglomeration Algebraic Multigrid Upscaling and Solvers) library, to produce multilevel preconditioners and solvers for H(curl) and H(div) formulations. ParELAG constructs hierarchies of compatible nested spaces, forming an exact de Rham sequence on each level. This allows the application of hybrid smoothers on all levels and AMS (Auxiliary-space Maxwell Solver) or ADS (Auxiliary-space Divergence Solver) on the coarsest levels, obtaining complete multigrid cycles. Numerical results are presented, showing the parallel performance of the proposed methods. As a part of the exposition, this paper demonstrates some of the capabilities of ParELAG and outlines some of the components and procedures within the library.

Section of the proceedings of FEniCS 2021

Recent years have seen a massive explosion of datasets across all areas of science, engineering, technology, medicine, and the social sciences. The central questions are: How do we optimally learn from data through the lens of models? And how do we do so taking into account uncertainty in both data and models?These questions can be mathematically framed as Bayesian inverse problems. While powerful and sophisticated approaches have been developed to tackle these problems, such methods are often challenging to implement and typically require first and second order derivatives that are not always available in existing computational models.We present an extensible software framework MUQ-hIPPYlib that overcomes this hurdle by providing unprecedented access to state-of-the-art algorithms for deterministic and Bayesian inverse problems. MUQ provides a spectrum of powerful Bayesian inversion models and algorithms, but expects forward models to come equipped with gradients/Hessians to permit large-scale solution. hIPPYlib implements powerful large-scale gradient/Hessian-based solvers in an environment that can automatically generate needed derivatives, but it lacks full Bayesian capabilities. By integrating these two libraries, we created a robust, scalable, and efficient software framework that realizes the benefits of each to tackle complex large-scale Bayesian inverse problems across a broad spectrum of scientific and engineering areas.Our ultimate goal is to make these advanced inversion capabilities available to a broader scientific community, to provide an environment that accelerates scientific discovery.

Recent years have seen a massive explosion of datasets across all areas of science, engineering, technology, medicine, and the social sciences. The central questions are: How do we optimally learn from data through the lens of models? And how do we do so taking into account uncertainty in both data and models?These questions can be mathematically framed as Bayesian inverse problems. While powerful and sophisticated approaches have been developed to tackle these problems, such methods are often challenging to implement and typically require first and second order derivatives that are not always available in existing computational models.We present an extensible software framework MUQ-hIPPYlib that overcomes this hurdle by providing unprecedented access to state-of-the-art algorithms for deterministic and Bayesian inverse problems. MUQ provides a spectrum of powerful Bayesian inversion models and algorithms, but expects forward models to come equipped with gradients/Hessians to permit large-scale solution. hIPPYlib implements powerful large-scale gradient/Hessian-based solvers in an environment that can automatically generate needed derivatives, but it lacks full Bayesian capabilities. By integrating these two libraries, we created a robust, scalable, and efficient software framework that realizes the benefits of each to tackle complex large-scale Bayesian inverse problems across a broad spectrum of scientific and engineering areas.Our ultimate goal is to make these advanced inversion capabilities available to a broader scientific community, to provide an environment that accelerates scientific discovery.

Non-smooth regularization is widely used in image reconstruction to eliminate the noise while preserving subtle image structures. In this work, we investigate the use of proximal Newton (PN) method to solve an optimization problem with a smooth data-fidelity term and total variation (TV) regularization arising from image reconstruction applications. Specifically, we consider a nonlinear Poisson-modeled single-energy X-ray computed tomography reconstruction problem with the data-fidelity term given by the I-divergence. The PN algorithm is compared to state-of-the-art first-order proximal algorithms, such as the well-established fast iterative shrinkage and thresholding algorithm (FISTA), both in terms of the number of iterations and time to solutions. We discuss the key factors that influence the performance of PN, including the strength of regularization, the stopping criterion for both sub-problem and main-problem, and the use of exact or approximated Hessian operators.

Slide to accompany 1 minute lightning talk at the NSF SI2 PI meeting. Our project is "Integrating Data with Complex Predictive Models under Uncertainty: An Extensible Software Framework for Large-Scale Bayesian Inversion"The goal of our project is to enable the widespread use of modern Bayesian inference algorithms by developing a software framework and a sustainable community of users and developers. Our framework is built on top of two major software packages: the MIT Uncertainty Quantification Library (MUQ), which provides tools for advanced statistical modeling and sampling approaches, as well as hIPPylib, which provides state-of-the-art tools for large-scale deterministic inverse problems built from partial differential equations.

This work describes a domain embedding technique between two non-matching meshes used for generating realizations of spatially correlated random fields with applications to large-scale sampling-based uncertainty quantification. The goal is to apply the multilevel Monte Carlo (MLMC) method for the quantification of output uncertainties of PDEs with random input coefficients on general, unstructured computational domains. We propose a highly scalable, hierarchical sampling method to generate realizations of a Gaussian random field on a given unstructured mesh by solving a reaction-diffusion PDE with a stochastic right-hand side. The stochastic PDE is discretized using the mixed finite element method on an embedded domain with a structured mesh, and then the solution is projected onto the unstructured mesh. This work describes implementation details on how to efficiently transfer data from the structured and unstructured meshes at coarse levels, assuming this can be done efficiently on the finest level. We investigate the efficiency and parallel scalability of the technique for the scalable generation of Gaussian random fields in three dimensions. An application of the MLMC method is presented for quantifying uncertainties of subsurface flow problems. We demonstrate the scalability of the sampling method with non-matching mesh embedding, coupled with a parallel forward model problem solver, for large-scale 3D MLMC simulations with up to $1.9\cdot 10^9$ unknowns., To appear in Numerical Linear Algebra with Applications

2017 NSF SI2 PI meetingFeb. 21-22, Arlington, VANSF SI2-SSI project:Integrating Data with Complex Predictive Models under Uncertainty: An Extensible Software Framework for Large-Scale Bayesian Inversion

Numerical solution of the system of partial differential equations arising from the Reynolds-Averaged Navier-Stokes (RANS) equations with $k-\epsilon$ turbulence model presents several challenges due to the advection dominated nature of the problem and the presence of highly nonlinear reaction terms. State-of-the-art software for the numerical solution of the RANS equations address these challenges by introducing non-differentiable perturbations in the model to ensure numerical stability. However, this approach leads to difficulties in the formulation of the higher-order forward/adjoint problems, which are needed for scalable Hessian-based uncertainty quantification (UQ) methods. In this note, we present the construction of a UQ-ready flow solver, i.e., one that is smoothly differentiable and provides not only forward solver capabilities but also adjoint and higher derivatives capabilities.

Current state of the art preconditioners for the reduced Hessian and the Karush-Kuhn-Tucker (KKT) operator for large scale inverse problems are typically based on approximating the reduced Hessian with the regularization operator. However, the quality of this approximation degrades with increasingly informative observations or data. Thus the best case scenario from a scientific standpoint (fully informative data) is the worse case scenario from a computational perspective. In this paper we present an augmented Lagrangian-type preconditioner based on a block diagonal approximation of the augmented upper left block of the KKT operator. The preconditioner requires solvers for two linear subproblems that arise in the augmented KKT operator, which we expect to be much easier to precondition than the reduced Hessian. Analysis of the spectrum of the preconditioned KKT operator indicates that the preconditioner is effective when the regularization is chosen appropriately. In particular, it is effective when the regularization does not over-penalize highly informed parameter modes and does not under-penalize uninformed modes. Finally, we present a numerical study for a large data/low noise Poisson source inversion problem, demonstrating the effectiveness of the preconditioner. In this example, three MINRES iterations on the KKT system with our preconditioner results in a reconstruction with better accuracy than 50 iterations of CG on the reduced Hessian system with regularization preconditioning., Comment: 30 pages, 4 figures, 1 table. Accepted for publication in SIAM Journal on Scientific Computing (SISC)